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SECTION  I 


Purpose  and  Scope 

The  purpose  of  the  spall  research  described  herein  is  defined  in  the 
Statement  of  Work  for  Air  Force  Office  of  Scientific  Research  Contract 
Number  F49620-81-C-0066,  viz: 

1.  Perform  Analytical  Studies  and  Finite  Difference  Calculations. 
Formulate  and  perform  one-  and  two-dimensional  analytical  studies  and 
calculations  to  obtain  insight  into  spall  mechanisms  and  important 
material  parameters.  Include  variations  in  wave  type  (e.g.  reflected  P 
and  S  waves,  horizontally  traveling  P  and  S  waves,  etc.),  tensile  failure 
critera  (e.g.,  maximum  principal  tension,  hydrostatic  tension,  extensional 
strain)  and  tensile  strength.  Include  variations  in  material  strength  if 
this  parameter  appears  to  be  important  in  spall  phenomena. 

2.  Develop  a  Spall  Prediction  Technique.  Using  the  results  of 

Task  1,  plus  previous  field  data  analysis,  develop  a  method  for  predicting 
spall  characteristics  in  both  single  and  multiburst  near-surface 
experiments.  The  prediction  technique  shall  include  the  prediction  of 
spall  radius,  depth,  initial  conditions,  and  post-spall  behavior,  all  as  a 
function  of  yield,  height  of  burst,  geology  and  material  properties.  The 
extension  of  the  method  for  prediction  of  spall  for  high  yield  nuclear 
bursts  shall  be  included. 

As  the  research  progressed,  two  things  became  clear: 

1.  Soil  spall  data  from  both  chemical  and  nuclear  explosive  tests  is 
limited  in  quantity  and  quality;  and 


APPLIED  RE/EAPCH  A//OCIATEf.mC. 


2.  A  fundamental  understanding  of  soil  spall  requires  an  effective 
stress  approach,  which  explicitly  treats  both  the  soil  skeleton  stress- 
strain  response  and  the  transient  effect  of  flowing  pore  fluid. 

Consequently,  relatively  more  attention  was  given  to  soil  tensile 
stress-strain  behavior,  and  also  particularly  to  the  pore  air  effect  than 
originally  planned,  and  somewhat  less  effort  was  spent  on  numerical 
calculations  (because  satisfactory  material  models  were  not  available). 

The  empirical  analysis  of  single  burst,  near-surface  explosion-induced 
soil  and  rock  spall  was  exhaustive,  and  led  to  the  single-burst  prediction 
technique  presented  herein.  However,  the  volume  of  test  data  for  multiple 
burst  spall  is  so  small,  and  the  material  behavior  influence  so  strong, 
that  a  multiple  burst  spall  prediction  technique  requires  a  combined 
empirical  and  theoretical  approach,  and  must  await  further  research. 

Section  III  of  this  report  describes  the  theoretical  basis  and  results 
of  one-dimensional  finite  difference  calculations  of  vertical  wave 
propagation  in  a  tensilely  weak,  hysteretic  material  with  a  free  surface, 
using  the  computer  program  STEALTH  ID.  Multiple  spall  features  are 
predicted  by  these  total  stress  calculations,  which  demonstrate  the 
ability  of  STEALTH  ID  to  both  include  gravity  in  propagation  problems,  and 
to  predict  spall  in  tensilely  weak  material  with  given  stress-strain 
behavior. 

Section  IV  reviews  previous  tests  and  analyses  of  the  pore  air  effect, 
and  presents  several  new  analyses  which  show  that  pore  air  expansion  and 
pore  air  flow  are  not  separate  spall  mechanisms,  but  rather  two  related 
aspects  of  the  same  mechanism.  Nonlinear  differential  equations  are 
derived  in  a  series  of  appendices,  for  both  isothermal  and  adiabatic 
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transient  pore  air  flow  through  a  rigid,  porous,  isotropic  solid.  Finite 
difference  approximations  to  the  above  differential  equations  are  also 
developed,  some  in  a  pseudo-linear  form  simple  enough  for  even  hand 
computation.  Two  simple,  discrete-element  models  of  the  pore  air  effect 
are  constructed  and  analyzed,  one  without  and  the  other  with  mass. 

Because  these  models  are  linear,  their  response  can  be  represented  in 
closed  form,  which  permits  parametric  studies.  Both  models  exhibit  rapid, 
undrained  response  (compression)  to  a  suddenly  applied  constant  surface 
fluid  pressure,  followed  by  gradual  recovery  (expansion)  of  the  soil 
skeleton  as  additional  pressurized  fluid  flows  into  the  pore  spaces.  The 
model  with  mass  is  a  fourth  order  model,  in  the  sense  that  its  closed  form 
analysis  requires  the  solution  of  a  fourth  order,  linear,  ordinary 
differential  equation  with  constant  coefficients.  This  involves  the 
solution  of  a  fourth  order  algebraic  equation,  with  one  root  several 
orders  of  magnitude  larger  than  the  rest.  Although  the  closed  form 
solution  of  this  equation  was  tedious,  it  was  pursued  because  it  provides 
considerable  insight  into  the  two  principal  modes  of  physical  behavior  of 
real  soil,  viz.  undrained  (composite)  behavior  and  drained  (separate) 
behavior  of  the  two  phases,  and  their  relative  rates  of  volumetric 
deformation.  Finally,  a  general  set  of  coupled,  nonlinear,  second  order 
partial  differential  equations  is  derived  which  simultaneously  describe 
both  the  propagation  and  diffusion  phenomena  which  comprise  the  pore  air 
effect.  These  equations  are  general  enough  to  accommodate  any  soil 
effective  stress-strain-strength  relationship.  Although  written  for  one¬ 
dimensional  motion,  they  can  readily  be  generalized  for  two-  or  three- 
dimensional  motion. 
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Sections  V,  VI  and  VII  comprise  an  integrated,  sequential  treatment  of 
explosion-induced  negative  (gage)  airblast  pressure  for  both  chemical  and 


nuclear  near-surface  bursts;  an  empirical  analysis  of  all  available  soil 
and  rock  spall  data  from  single,  near-surface  explosions,  both  chemical 
and  nuclear;  and  an  empirically-based  predictive  method  for  soil  and  rock 
spall  from  single,  near-surface  explosions,  both  chemical  and  nuclear. 
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SECTION  II 


Introduction 

The  study  of  near-surface  explosion  induced  soil  spall  is  basically  a 
study  of  dynamic  soil  stress-strain  behavior,  and  requires  explicit 
recognition  of  the  particulate,  multiphase  nature  of  soil.  This  does  not 
mean  that  total  stress  approaches  are  not  useful,  but  rather  that  their 
basis  and  interpretation  are  rooted  in  the  concept  of  effective  stress. 
Thus  the  finite  difference  calculations  in  Section  III,  and  the  empirical 
analyses  in  Sections  VI  and  VII  use  a  total  stress  approach,  but  the 
analyses  of  the  pore  air  effect  in  Section  IV  use  an  effective  stress 
approach,  and  help  explain  the  results  presented  in  the  other  sections. 

Much  work  remains  to  be  done  to  develop  a  reasonably  thorough 
understanding  of  soil  spall  and  its  effect  on  near-surface  explosion 
induced  ground  motions.  Much  of  that  work  is  of  a  calculational  nature, 
involving  parametric  studies  using  the  equations  developed  herein,  and 
comparing  the  calculated  results  with  the  results  of  predictions  made 
using  empirical  techniques  also  developed  herein.  From  this  point  on  the 
further  development  of  spall  calculation  and  empirical  prediction 
techniques  can  best  proceed  in  tandem. 
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SECTION  III 


Numerical  Studies  of  Spall  Under  Explosive  Loading 
A.  Purpose  of  the  Calculations 

Mathematical  simulation  and  prediction  of  explosive  events  in  soil 
require  adequate  treatment  of  tensile  (spall)  phenomena.  Spallation, 
whether  caused  by  layering  effects,  wave  interaction  with  the  ground 
surface,  or  some  other  phenomena,  will  significantly  affect  measured 
ground  motions  (especially  at  later  times)  and  therefore  should  be  modeled 
in  a  calculation.  In  order  to  illustrate  the  generation  of  spall  and 
examine  the  effect  of  varying  physical  parameters  such  as  tensile 
strength,  a  series  of  one-dimensional  finite  difference  calculations  was 
performed  for  this  study. 

Figure  3.1  shows  two  situations  which  may  lead  to  spall.  In 
Figure  3.1a,  shock  waves  produced  by  an  explosion  at  or  above  the  earth's 
surface  are  propagated  much  faster  in  an  underlying  layer  of  rock  than  in 
the  surface  layer  of  soil.  As  a  result,  a  compressive  headwave  moves 
upward  into  the  undisturbed  soil.  The  angle  from  the  vertical  (e)  at 
which  this  wave  propagates  may  be  estimated  as: 

e  =  arcsin  (Cs/Cr)  (3.1) 

where  Cs  and  Cr  are  the  confined  wavespeeds  in  soil  and  rock,  respectively. 
In  a  situation  where  the  lower  layer  has  a  much  higher  wavespeed  than  the 
upper  one  (i.e. ,  a  high  impedance  mismatch),  the  headwave  will  travel 
nearly  vertically.  Because  of  this,  a  one-dimensional,  uniaxial 
idealization  is  appropriate,  as  shown  in  Figure  3.1c.  Such  a  1-D 
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calculational  model  of  spall  may  also  be  adequate  for  the  case  of  a  buried 
burst,  as  shown  in  Figure  3.1b.  Here  the  primary  waves  move  vertically 
through  undisturbed  soil  and  reflect  off  the  ground  surface  as  tensile 
waves  which  cause  spall. 

The  one-dimensional  model  used  for  these  calculations  (Figure  3.1c) 
may  be  considered  to  be  an  adequate  simplification  for  studying  many  real 
cases  of  spall.  It  is  true,  however,  that  certain  occurences  of  spall  can 
be  accurately  simulated  only  in  two  or  three  dimensions.  Such  cases 
require  more  sophisticated  treatment  of  tensile  behavior  and  geometric 
factors. 

B.  Computer  Code  Description 

There  are  many  one-,  two-,  and  three-dimensional  finite  difference 
wave  propagation  codes  currently  available.  The  particular  code  chosen  to 
perform  this  set  of  calculations  was  STEALTH  [Hofmann  (1978)].  STEALTH 
1-D,  2-D,  and  3-D  are  a  set  of  "user-oriented"  codes,  and  have  been 
applied  to  many  different  kinds  of  dynamic  problems. 

STEALTH  is  attractive  for  studing  spall  because  of  its  ability  to 
accept  different  constitutive  relationships  and  to  allow  the  user  to 
impose  a  wide  variety  of  boundary  conditions.  The  generality  of  STEALTH, 
although  it  requires  a  very  large  code,  is  also  one  of  its  strong  points. 
For  example,  a  purely  mechanical  spall  model  worked  in  one-dimension  could 
be  broadened  to  include  pore  fluid  effects  and  still  be  implemented  in 
STEALTH  1-D  or  2-D,  sparing  any  further  development  necessitated  by  code 
differences. 

C.  Description  of  Calculational  Set-Up 

The  calculations  for  this  study  were  performed  under  uniaxial  strain 
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conditions.  A  27  meter  grid,  consisting  of  160  zones  was  used,  as  shown 
in  Figure  3.2.  The  idealized  site  has  a  layered  geology  with  12  m  of 
"soil"  over  15  m  of  "rock".  Actually,  the  soil  is  simply  a  bilinear, 
hysteretic  material  with  an  unloading  wave  speed  (1122  m/s)  equal  to  twice 
the  loading  wave  speed  (561  m/s).  The  rock  is  modeled  elastically  with 
equal  loading  and  unloading  speeds  (2960  m/s).  A  Pm^n  (minumum  bulk 
stress)  criterion  was  used  to  establish  the  soil  spall-threshold.  Because 
of  the  one-dimensional  nature  of  the  calculations,  this  is  equivalent  to  a 
vertical  principal  stress  cutoff  criterion.  The  tensile  cutoff  (P  ^n) 
was  varied  for  the  soil,  but  the  rock  was  not  allowed  to  fail  in  tension. 
Neither  material  was  allowed  to  fail  in  shear. 

Gravity  is  an  important  part  of  any  spall  situation  because  it  creates 
the  lg  dwell  signatures  which  are  a  primary  identifier  of  spall  in  data 
records.  In-situ  gravity  stresses  were  established  in  the  STEALTH  grid  by 
applying  a  body  force  equivalent  to  lg  over  the  entire  grid  with  a  rigid 
boundary  at  the  bottom.  The  resultant  oscillations  were  then  numerically 
damped  (in  STEALTH,  this  is  known  as  "dynamic  relaxation").  Within  100  ms 
a  static  condition  was  achieved.  Figure  3.3  shows  grid  velocities  at 
various  depths  during  this  relaxation  phase.  All  information  concerning 
the  grid  was  saved  at  200  ms  and  subsequent  calculations  were  restarted 
from  this  time.  Figure  3.4  shows  the  in-situ  geostatic  stresses  in  the 
grid  after  gravity  has  been  established.  Note  the  discontinuity  in  bulk 
stress  at  the  soil-rock  interface.  This  is  a  result  of  unequal  Kq 
values  for  the  two  materials.  This  may  or  may  not  be  indicative  of  actual 
in-situ  conditions  at  a  soil-rock  interface,  because  of  uncertain  geologic 
history. 
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The  boundary  condition  at  the  top  of  the  grid  for  ail  problems  was  a 
free  surface,  i.e. ,  zero  stress  (atmospheric  pressure  was  neglected).  In 
order  to  simulate  the  effect  of  an  infinite  rock  layer  below  the  soil,  a 
transmitting  (or  nonreflecting)  boundary  was  used  at  the  bottom  of  the 
grid  after  gravity  had  been  established.  A  transmitting  boundary  allowing 
incident  waves,  which  has  been  used  by  [Moriwaki  et.  al.  (1981)]  for  shear 
waves,  was  adapted  for  compressional  waves  and  the  effect  of  gravity. 
Basically,  the  velocity  of  the  boundary  point  is  calculated  based  on  a 
a  =  pcv  relationship.  Referring  to  Figure  3.5,  Newton's  Law  for  the 
boundary  element  may  be  expressed  as: 

m  x  =  EFg  (3.2) 

where 


=  boundary  mass  =  j  pax 


(3.3) 


EFg,  the  sum  of  the  live  forces  acting  on  the  boundary  element,  includes 
(o  -  og),  the  live  stress  in  the  zone  above  the  boundary,  and  a g,  the 
boundary  stress,  where 

aB  =  Pc( x  -  2X)  (3.4) 

The  quantity  c  is  the  compressional  wave  speed  of  the  imaginary 
infinite  media  below  the  boundary,  and  X  is  the  applied  incident  velocity 
at  the  boundary.  With  the  timing  of  the  calculation  as  shown  in 
Figure  3.5.,  and  using  superscripts  n  and  1  to  denote  new  and  old  times, 
respectively,  the  current  acceleration  and  velocities  are: 


(3.5b) 
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(3.5c) 


Substituting  in  Equation  (3.2)  for  m,  x,  and  ZFb  yields 


-1  ...wxn  -  x\  /  v  _  f,xn  +  x1, 


(7  pax) ( — 


-)  =  (a  -  ar)  -  pC  ( 


-  (X°  +  X1)] 


(3.6a) 


and  solving  for  the  desired  boundary  velocity,  xn  yields 


(*"  -  *’)  -  -  s&i"  *  x')  -  2£ii(x"  ♦  x') 


(3.6b) 


cat,  +  2At(g  -  ctg)  _  2cat  (vT)  +  yh 

AXJ _ PAX _ AX  *  *  J 

(1  +  ci|) 


(3.7) 


Using  the  boundary  condition  defined  by  Equation  3.7,  it  is  possible 
to  propagate  a  pulse  up  into  undisturbed  material  while  simultaneously 
allowing  reflected  waves  to  pass  down  through  the  bottom  of  the  grid  and 
out  of  the  problem.  Figure  3.6  shows  characteristic  planes  for  two  cases, 
one  with  a  pressure  boundary  at  the  bottom  and  another  with  the 
transmitting  boundary.  It  is  apparent  that  the  latter  case.  Figure  3.6b, 
with  no  increase  in  grid  size,  more  closely  models  actual  conditions,  such 
as  those  shown  in  Figure  3.1.  Therefore  the  transmitting  boundary  allows 
analysis  of  spall  motions,  unobstructed  by  peripheral  reflections. 

The  velocity  pulse  used  to  generate  spall  in  the  top  soil  layer  is 
shown  in  Figure  3.2.  It  consists  of  a  linear  rise  with  subsequent 
sinusoidal  decay  to  zero.  The  peak  velocity,  Vp^,  was  varied,  while 
maintaining  the  pulse  duration  at  20  ms. 

D.  Discussion  of  Results 

After  establisning  suitable  boundary  conditions,  several  calculations 
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were  performed  to  meet  the  following  objectives: 

i)  Examine  the  case  of  no  spall,  where  the  soil  material  is  not 
allowed  to  separate  in  tension; 

ii)  Examine  a  baseline  spall  problem,  for  which  the  general 
characteristics  of  spall  can  be  examined  and  against  which  the 
results  of  parameter  variations  can  be  compared; 

iii)  Assess  the  effect  of  load  pulse  magnitude  on  spall; 

iv)  Assess  the  effect  of  soil  tensile  strength  on  spall; 

v)  Determine  the  importance  of  numerical  noise  in  the  calculations 
which  may  produce  spall  motions;  and 

vi)  Examine  a  case  in  which  spall  blocks  are  forced  back  together 
("forced  rejoin"),  in  order  to  ascertain  the  basic 
characteristics  of  the  resulting  waveforms. 

Table  3.1  is  a  matrix  of  calculations  with  values  for  the  parameters 
which  were  varied,  viz;  tensile  cutoff  and  peak  incident  velocity  at  the 
bottom  of  the  grid.  In  addition,  for  all  the  calculations  except  No.  6, 
the  tensile  stress  in  a  zone  was  required  to  exceed  PMIN  for  three 
cycles  before  spall  was  allowed  to  occur.  All  problems  were  started  with 
an  established  gravity  stress  field  at  time  =  0  and  were  run  for  250  ms. 

The  results  from  Calculation  No.  1  illustrate  what  the  waveforms  in  a 
material  with  a  large  tensile  strength  would  look  like  for  the  loading  of 
interest.  The  material  remains  intact  and  the  oscillations  can  be  easily 
predicted.  Figure  3.7  shows  the  calculated  stress  and  velocity  time 
histories  at  a  depth  of  0.5  m.  (Note  the  sign  conventions,  which  will  be 
used  throughout  this  section.)  Below  each  waveform  the  characteristic 
plane  for  the  problem  is  shown.  As  reflected  energy  passes  through  the 
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rock  and  out  of  the  grid,  the  motions  eventually  damp  out  to  zero. 

When  the  tensile  strength  of  the  soil  is  reduced  to  a  reasonable 
level,  motions  are  dramatically  altered.  Figure  3.8  shows  a  comparison 
between  the  cases  of  spall  (PMjN  =  0.10  MPa)  and  no  spall  (PMIN  = 

100  MPa). 

There  are  several  spall-significant  events  in  the  calculation  which 
contribute  to  the  waveform  differences  seen  in  Figure  3.8.  The  following 
is  a  summary  of  these  events  in  chronological  order,  along  with  a  list  of 
the  more  pertinent  parameters  controlling  them: 

i)  The  upward  propagating  compressive  pulse  hits  the  boundary 
between  soil  and  rock  (time  =  5  ms).  Part  of  the  wave  reflects 
back  into  the  rock  as  a  tensile  wave,  and  part  transmits  through 
as  a  compressive  wave.  At  this  point,  the  soil  separates  from 
the  rock  (see  Figure  3.9). 

parameters:  PMJN  (soil  and  rock) 
cu  (soil  and  rock) 

ii)  The  initial  compressive  wave  travels  upward,  compressing  the 
hysteretic  soil  material  (see  Figure  3.10).  Note  that  because  a 
net  upward  displacement  is  imposed  on  the  bottom  of  the  rock 
layer  by  the  incident  velocity  pulse,  the  final  position  of  a 
soil  node  is  not  dependent  only  on  soil  compression.  When  the 
free  surface  is  encountered  (at  approximately  26  ms),  the  wave 
reflect;  as  a  tensile  wave.  At  the  depth  at  which 

Reflected  <  PMIN 

spall  first  occurs. 

parameters:  tsl  =  time  of  first  spall 
d$j  =  depth  of  first  spall 
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iii)  The  tensile  wave  continues  to  propagate  downward  at  the  soil 


unloading  wavespeed,  and  may  or  may  not  create  more  spall  planes 
depending  on  the  material  model  being  used.  Eventually  it 
reaches  the  free  surface  which  has  been  created  between  the  rock 
and  soil.  At  this  time  {approximately  37  ms)  the  entire  soil 
block  realizes  it  is  undergoing  projectile  motion.  The  tensile 
pulse  must  reflect  as  a  compression  wave. 

parameters:  n  =  number  of  initial  spall  planes  formed 

tc  ,  =  times  of  subsequent  initial  spall 

plane  formation 
d  ~  =  depths  of  above 

b}  u  •  •  •  1 1 

iv)  There  is  a  period  of  motion  with  constant  acceleration 
(projectile  motion,  lg  dwell,  freefall,  etc.)  for  the  soil. 

Energy  trapped  within  the  discrete  spall  blocks  moves  between  the 
free  surfaces  of  the  blocks.  This  causes  small  oscillations 
during  the  lg  dwell,  whose  frequency  depends  directly  on  spall 
block  size  (see  Figure  3.8).  The  magnitude  of  these  trapped 
waves  may  be  high  enough  in  some  blocks  to  cause  further  spall. 

parameters:  t^  =  duration  of  lg  dwell 
(varies  with  depth) 

v)  Rejoin  of  spalled  material  begins  at  the  rock  surface.  The  time 
at  which  rejoin  occurs  for  each  node  is  roughly  determinable  from 
the  initial  velocity  (xQ)  and  displacement  (xQ)  at  the  time 
spall  occurs.  The  equation  of  motion  for  constant  acceleration 
(x  =  g)  is: 

x  =  % +  v +  ¥t2  (3’8) 
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where  t  =  time  from  spall.  Solving  for  the  duration  of  spall, 
td,  yields 


In  Equation  (3.9),  xd  is  the  particle  position  at  the  time  of 
rejoin,  and  depends  on  the  compressibi 1 ity  of  the  soil  and  the 
incident  velocity  pulse.  For  estimating  dwell  time,  xrf  may  be 
set  to  zero,  i.e.,  the  pre-spall  position.  Figure  3.11 
illustrates  the  calculation  of  dwell  time  at  the  5  m  depth, 
parameters:  tr  =  t  +  td  =  time  of  rejoin 

vi)  Rejoin  progresses  from  the  bottom  of  the  soil  upwards.  When  two 
spall  blocks  impact,  a  compressive  signal  travels  both  upwards 
and  downwards.  An  upward  pulse  encounters  a  free  surface 
(because  the  block  above  has  not  rejoined  yet)  and  reflects. 

This  reflection  may  cause  further  spall  in  the  soil.  The 
downward  pulse  reflects  and  transmits  at  the  rock  boundary.  A 
negative  pulse  (due  to  the  impedance  mismatch)  follows  this 
reflection  up  and  may  also  cause  further  spall.  As  a  result, 
spallation  of  zones  accompanies  rejoin  and  the  number  of  spall 
planes  increases  due  to  impact,  creating  somewhat  of  a 
"shattering"  effect.  Rejoin  propagates  slower  as  it  approaches 
the  surface  because  the  upper  soil  zones  have  been  lofted 
somewhat  higher. 

parameters:  vr  -  rejoin  (or  impact)  velocities 

vi i)  Subsequent  spall  motion  occurs  and  may  be  referred  to  as 
secondary  or  multiple  spall.  During  this  phase,  spall  planes 
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are  simultaneously  opening  and  closing.  Figure  3.12  shows  the 
formation  of  several  spall  planes  due  to  the  first  reflection, 
and  subsequent  spall  plane  formation  during  the  rejoin  and 
secondary  spall  phases.  As  rejoin  energy  is  transmitted  through 
the  soil/rock  boundary  and  through  the  nonreflecting  boundary, 
motions  damp-out  and  the  grid  elements  come  to  rest, 
parameters:  t  =  times  of  multiple  spall 

t  =  times  of  multiple  spall  rejoin 
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SECTION  IV 


The  Pore  Air  Effect 


A.  Introduction 

Spall  in  dynamically  loaded  granular  soil  occurs  when  soil  masses 
lose  contact  with  each  other  and  undergo  ballistic  motion  (free  fall). 

An  important  cause  of  spall  near  the  ground  surface,  in  the  presence  of 
explosion-induced  airblast  pressure  is  the  pore  air  effect.  Although 
sometimes  discussed  as  a  parallel  phenomenon  with  spall,  the  pore  air 
effect  is  best  considered  as  a  spall  mechanism  [Merkle  (1980:24)]. 

The  pore  air  effect  is  a  consequence  of  soil  being  both  particulate 
and  multiphase  [Lambe  and  Whitman  (1969:18-19)].  Under  the  action  of  a 
pore  pressure  gradient,  soil  pore  fluid  flows  in  the  direction  opposite 
to  that  of  the  pore  pressure  gradient.  An  increase  in  the  mass  of  pore 
fluid  stored  in  a  saturated  soil  element  occurs  if  the  soil  void  volume 
increases,  and/or  the  pore  fluid  mass  density  increases.  In  general, 
both  can  happen.  If  inertia  is  neglected,  and  the  equations  for  pore 
fluid  flow  and  storage  are  combined,  the  result  is  a  diffusion  equation 
which  can  be  written  with  pore  pressure  as  the  unknown. 

If  the  pore  fluid  flow  equation  is  assumed  to  be  linear,  the  pore 
fluid  incompressible,  and  the  soil  skeleton  linearly  elastic,  the  result 
is  Terzaghi's  consolidation  equation  [Terzaghi  (1943:265);  Taylor 
(1948:^25) ]. 

If  the  pore  fluid  flow  equation  is  again  assumed  to  be  linear,  the 
pore  fluid  (air)  to  undergo  isothermal  compression,  and  the  soil  skeleton 
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to  be  rigid,  the  result  is  a  nonlinear  diffusion  equation  [Zernow  et  al 
(1973:13);  Muskat  (1937);  Carman  (1956)]. 

If  the  soil  particles  are  assumed  to  be  in  a  spalled  condition,  i.e., 
not  in  contact,  and  the  pore  air  to  undergo  adiabatic  compression,  the 
resulting  some  velocity  (which  applies  to  a  propagation,  rather  than  a 
diffusion  problem)  is  that  for  a  heavy  gas  [Ullrich  (1978:19);  Merkle 
(1980:35)].  In  this  model  the  mass  of  the  heavy  gas  is  supplied  entirely 
by  soil  narticles,  and  the  compressibi 1 ity  entirely  by  pore  air.  There 
is  no  viscous,  diffusive  flow  in  this  model,  but  there  is  compressible 
flow  in  the  sense  of  compressive  wave  particle  motion. 

The  above  pore  air  models  and  other  more  general  models  are  discussed 
in  detail  below.  Emphasis  is  on  developing  equations  for  pore  air 
pressure  prior  to  spall,  to  gain  insight  into  how  pore  air  pressure  can 
cause  spall. 

Two  pore  air  phenomena  have  been  proposed  as  possible  causes  of  spall: 
pore  air  expansion,  and 
pore  air  flow. 

Saying  that  pore  air  expansion  causes  spall  can  be  misleading,  because 
adiabatic  pore  air  expansion  is  accompanied  by  a  pore  pressure  decrease, 
which,  by  itself,  would  cause  an  increase  in  effective  hydrostatic  stress 
carried  by  the  soil  skeleton,  resulting  in  compression  rather  than 
expansion  of  the  soil  skeleton.  What  actually  happens  when  airblast 
pressure  acts  on  a  soil  surface  is  that  air  flows  into  the  soil, 
compressing  the  pore  air  near  the  loaded  surface.  The  compressed  pore 
air  has  a  tendency  to  expand  rapidly,  especially  when  the  surface 
airblast  pressure  decreases,  and  this  rapid  expansion  produces  a  flow  of 
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pore  air  toward  the  soil  surface.  The  pore  air  flow  in  turn  causes  a 
seepage  or  drag  force  to  act  on  soil  particles  in  the  flow  regime,  which 
tends  to  reduce  the  effective  stresses  produced  by  gravity  and  eventually 
may  cause  spall.  Thus  it  is  pore  air  expansion  and  flow,  together,  which 
can  cause  spall.  A  profitable  way  to  view  the  pore  air  effect  is  to 
consider  the  pore  air  seepage  force  (the  negative  of  the  pore  pressure 
gradient)  as  a  body  force  acting  on  the  soil  skeleton,  the  deformation  of 
which  is  governed  by  effective  stress.  The  effective  stress  carried  by 
the  soil  skeleton,  and  the  pore  air  pressure  are  always  such  that  the 
soil  void  volume  and  the  pore  air  volume  are  equal. 

B.  Previous  Investigations 

1 .  Hampton  (1964) 

Hampton  performed  a  shock  tube  study  of  dynamic  pore  air  pressure  in 
three  dry  soils  subjected  to  surface  airblast  pressure.  The  loaded  soil 
surface  was  confined  by  a  screen  (presumably  to  prevent  lofting  or 
spall),  because  interest  centered  on  the  rate  of  attenuation  of  peak  pore 
air  pressure  in  soil,  as  it  affects  the  design  of  buried  structures, 
foundations  for  surface  and  buried  structures,  and  model  experiments  on 
soil-structure  interaction.  In  particular,  Hampton  was  concerned  with 
the  depth  to  which  pore  air  pressure  induced  in  soil  will  penetrate,  the 
distribution  of  pore  air  pressure  with  depth,  and  the  relative  velocity 
of  pore  air  pressure  propagation  and  effective  stress  waves.  The  study 
was  entirely  experimental,  but  with  a  final  recommendation  that  an 
analytical  theory  of  pore-air  pressure  propagation  in  soil  subjected  to 
an  air  shock  wave  be  developed.  Hampton  measured  the  physical 
permeability,  k  ,  of  the  three  soils  tested,  with  the  following  average 
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results: 


pea  gravel:  2.42  X  10~6  CM2 
20-30  Ottawa  sand:  6.11  X  10~7  CM2 
silty  sand:  9.81  X  10"10  CM2 

Physical  permeability  is  used  in  the  Darcy  flow  equation  in  the  form 


where 

v  =  flow  velocity,  CM/SEC 

.  p 

kQ  =  physical  permeability,  CM 

u  =  viscosity,  DYNE  SEC/CM2 

P  =  pressure,  DYNE/CM2 

s  =  distance,  CM 

The  value  of  M  for  air  is  [Vennard  (1954:7)] 


yAIR 


3.77  X  10~7 

FT^ 


(3.77  X  10  7)(4. 448222  X  105)  =  1>81  x  10~4  DYNE  SEC 
(30.48)2  ’  CM2 


Equation  (4.1)  is  often  written  in  the  form 


v  =  -B 


dP 
1  ds 


(4.2) 


where  the  effective  permeability  coefficient,  B-^,  is  given  by  the 
expression 


For  air  flowing  through  standard  20-30  Ottawa  sand,  the  above  values  yield 
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B1  - 


6.11  X  10 


-7 


1.81  X  10 


^r 


=  3.38  X  10 


-3 


CNf 

DYNE  S£C 


2.  Zernow  et  al  (1973) 

Zernow  et  al  also  performed  a  shock  tube  study  of  Number  30  silica 
sand  subjected  to  surface  airblast  overpressure.  However,  their 
principal  objectives  were  to  determine  the  amount  of  near-surface  soil 
lofted  (or  spalled)  by  the  airblast-induced  "reverse  percolation" 
process,  and  to  measure  the  resulting  soil  particle  motions.  Their  study 
was  motivated  by  the  possibility  of  dust  cloud  formation  by  a 
near-surface  nuclear  detonation,  caused  by  a  combination  of  lofting  and 
airblast  sweep-up.  The  possibility  of  such  a  phenomenon  had  been 
suggested  by  Brode  in  1971. 

A  theoretical  description  of  the  lofting  process  was  also 
accomplished,  based  on  a  linear  form  of  the  equation  for  isothermal 
diffusion  of  air  through  a  rigid,  porous  medium.  The  relevant  equations, 
which  were  not  derived  in  the  above  report,  are  derived  in  Appendix  A. 

The  basic  linear  diffusion  equation  employed  to  calculate  pore  air 
pressure  was 


3P  n9  P 


where  the  diffusion  coefficient  for  isothermal  flow  is 


(A.6) 


D  = 


V 


(A. 7) 


Although  the  shock  tube  soil  samples  were  of  finite  length  (35  IN  for  the 
4  IN  shock  tube  [Zernow,  Figure  9],  and  71  IN  for  the  8  IN  shock  tube 
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[Zernow,  Figure  22]),  the  boundary  and  initial  conditions  used  to  solve 
Equation  (A. 6)  were  those  for  a  half space  [Carslaw  and  Oaeger  (1959:64)]: 


P(o.t)  =  Poe_at 

(t>0) 

(4.4) 

P(»,t)  =  o 

(t>0) 

(4.5) 

P(x,0)  =  0 

(x>0) 

(4.6) 

The  rationale  for  Equation  (4.5)  was  that  the  standard  soil  samples  were 
supported  by  a  vented  end  plate  having  a  porosity  of  about  19  percent  for 
the  4  IN  shock  tube  and  16  percent  for  the  8  IN  shock  tube,  whereas  the 
sample  porosity  was  about  33  percent.  This  the  authors  argued  would 
simulate  the  flow  resistance  of  the  imaginary  soil  beneath  the  finite 
sample  [Zernow,  pp.  39  and  71].  Unfortunately,  two  facts  contradict  the 
above  argument: 

a)  No  attempt  was  made  to  simulate  the  pore  air  storage 
characteristics  of  the  imaginary  soil  beneath  the  finite  sample. 

b)  With  a  vented  end  plate  it  is  difficult  to  conclude  that  the  pore 
air  pressure  at  the  sample  bottom  could  have  been  anything  but 
close  to  zero  (gage) . 

Therefore  transient  solutions  to  Equation  (A. 6)  have  been  obtained  in 
Appendices  B  and  C  for  the  following  boundary  and  initial  conditions: 


Appendix  B  (vented  sample) 


P(0,t)  =  P0e'at 

(t>0) 

(4.4) 

P(l,t)  =  0 

( t>o) 

(B.2) 

P(x,0)  =  0 

( 0<  x<  1 ) 

(B.3) 
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(4.4) 


<• 


Appendix  C  (unvented  sample) 

P(0,t)  =  P0e"at  (t>0) 

-§(l,t)-0  (t>0)  (C.2) 

P  ( x  ,0 )  =  0  (Ckx<  1 )  (C.3) 

The  reason  for  deriving  the  above  solutions  was  to  give  the  linear 
diffusion  theory  as  realistic  a  test  as  possible.  If  the  linear  theory 
is  judged  inadequate  for  predicting  dynamic  pore  air  pressures,  it  should 
be  because  the  phenomenon  really  is  not  linear,  rather  than  because 
inappropriate  boundary  conditions  were  used  for  the  linear  theory.  The 
linear  theory  should  not  be  hastily  abandoned,  because  as  Zernow  et  al 
point  out,  it  serves  as  a  source  of  useful  physical  insight  in  terms  of 
parametric  variations. 

Figures  4.1,  4.2  and  4.3  show  dynamic  pore  air  pressure  isochrones 
for  a  vented  sample  of  finite  length,  subjected  to  a  decaying  exponential 
airblast  pressure  on  one  face.  Measured  peak  pore  air  pressures  and 
their  time  of  occurrence  are  shown  in  Table  4.1,  and  the  appropriate 
linear  diffusion  parameters  are  calculated  in  Table  4.2.  Calculated 
values  for  comparison  with  the  measured  values  shown  in  Table  4.1  are 
shown  in  Figures  4.4,  4.5  and  4.6  and  in  Table  4.3.  The  agreement 
between  measured  values  of  P^x  and  t ^  shown  in  Table  4.1,  and 
corresponding  calculated  values  shown  in  Table  4.3  is  not  good.  Because 
of  this  lack  of  agreement,  and  because  Zernow' s  measured  soil 
permeability  values  have  already  been  questioned  by  [Morrison  (1979a: 5) J, 
the  above  calculations  were  rerun  for  soil  Groups  I,  II  and  III  using  the 
effective  permeability  for  Ottawa  sand  calculated  following  Equation 
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(4.3)  (3.4  X  10-^  CM^/DYNE  SEC).  Since  this  value  is  one-fifth  that 
shown  in  Table  4.1  for  Groups  I,  II  and  III,  the  adjusted  value  for 
Group  IV  in  Table  4.4  was  also  taken  as  one-fifth  that  shown  in  Table  4.1 
for  Group  IV.  The  comparison  between  calculated  values  of  and 
tMAX  shown  in  Table  4.5  and  corresponding  values  shown  in  Table  4.1  is 
much  improved.  The  degree  of  agreement  could  be  further  improved  by 
decreasing  the  assumed  permeability  even  more  to  account  for  the  greater 
fraction  of  fine  sand  and  silt  size  particles  in  Zernow's  material  than 
in  standard  Ottawa  sand.  It  thus  appears  that  the  early  phase  of  the 
pore  air  effect  prior  to  spall  can  be  described,  at  least  approximately, 
by  a  linear  diffusion  model. 

Note  the  upward  pore  air  flow  near  the  surface  in  Figures  4.8d,  e 
and  f.  This  is  the  flow  which  can  cause  spall  (lofting),  and  the  depth 
of  spall  can  be  estimated  as  the  depth  to  the  point  at  which  the 
isochrone  slope  is  zero. 

Figures  4.13  through  4.18  show  the  results  of  a  linear  diffusion 
analysis  of  transient  pore  air  flow  in  an  unvented  sample,  using 
Equation  (C.29).  The  effect  of  the  impervious  boundary  at  x  =  1  is 
dramatic,  and  the  tendency  toward  upward  air  flow  much 
greater  than  for  a  vented  sample.  The  impervious  boundary  could  be  a 
water  table,  a  rock  layer,  or  even  a  clay  layer. 

Because  Zernow  et  al  did  not  obtain  satisfactory  agreement  between 
measured  pore  air  pressures  and  motions  and  those  predicted  using  linear 
diffusion  theory,  they  developed  a  more  accurate,  nonlinear  diffusion 
theory.  Their  nonlinear  theory  considers  the  pore  air  flow  to  be 
adiabatic  rather  than  isothermal,  and  recognizes  the  spatial  variation  of 
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pore  air  density.  The  equations  are  developed  in  Appendix  D.  Comparison 
of  measured  and  predicted  peak  pore  air  pressure  and  time  to  peak 
pr»?ssure  using  the  nonlinear  theory  was  still  not  satisfactory,  but 
halfspace  boundary  conditions  were  again  used  instead  of  those  for  a 
finite  layer.  Because  of  the  unsatisfactory  comparison,  the  authors 
recommended  further  development  of  the  nonlinear  model  to  eliminate  the 
constraints  of  a  rigid  medium  with  constant  properties. 

The  authors  also  stated  that  there  was  no  information  on  how 
permeability  varies  with  porosity.  However,  this  subject  has  been  fairly 
thoroughly  studied  [Taylor  (1948:111);  Leonards  (1962:121);  Lambe  and 
Whitman  (1969:283);  Mitchell  (1976:346)]. 

3.  Ullrich  (1978) 

Ullrich  examined  MISER'S  BLUFF  single  and  multiple  burst  data  to 
determine  the  limits  of  superposition  as  a  means  of  predicting  multiple 
burst  groundshock  response  using  single  burst  groundshock  data.  He 
concluded  that  superposition  of  single  burst  groundshock  records  yielded 
good  predictions  of  multiple  burst  groundshock  motions  at  all  depths 
outside  the  multiple  charge  array,  and  for  depths  greater  than  10  feet 
inside  the  array.  At  depths  less  than  10  feet  inside  the  array 
superposition  failed,  and  he  suggested  four  mechanisms  which  caused 
superposition  to  fail: 

a)  airblast  pressure  enhancement 

b)  soil  rebound 

c)  soil  tensile  failure  [spall] 

d)  pore  air  expansion 

Ullrich  used  the  term  pore  air  expansion  to  describe  "the  expansion, 
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during  the  negative  overpressure  phases,  of  air  initially  entrained  in 
the  ground",  and  concluded  that  this  mechanism  was  the  dominant  cause  of 
superposition  failure  in  MISER'S  BLUFF. 

Ullrich  noted  that  the  single  burst  ground  motion  records  used  in  the 
attempted  prediction  of  a  near  surface  multiple  burst  ground  motion 
response  near  the  center  of  a  charge  array  did  not  show  spall,  but  the 
corresponding  multiple  burst  ground  motion  record  did  show  spall.  This 
situation  points  out  the  existence  of  two  classes  of  spall  prediction 
problems  in  a  superposition  context: 

a)  how  to  combine  component  records  which  themselves  show  spall,  and 

b)  how  to  combine  component  records  which  themselves  do  not  show 
spall,  but  which  in  combination  will  (in  reality)  cause  spall. 

Ullrich  realized  that  the  only  fundamental  approach  to  such  nonlinear 
problems  is  with  an  accurate,  nonlinear  material  model,  and  he  proposed 
one  for  already  spalled  soil.  The  model  is  that  of  a  heavy  gas, 
consisting  of  disconnected  soil  particles  suspended  in  air.  His  analysis 
used  the  term  "piston  model",  but  in  fact  the  model  can  be  used  for  two- 
or  three-dimensional  analyses  because  it  is  basically  a  volumetric 
model.  Ullrich's  equations  are  derived,  using  conventional  soil 
mechanics  nomenclature,  in  Appendix  E.  Although  Ullrich's  model  has  been 
called  a  "no  flow"  model,  it  does  yield  particle  motion  associated  with  a 
compressive  wave.  It  is  a  propagation,  not  a  diffusion  model,  and  does 
not  treat  viscous  flow  of  air  through  the  soil  skeleton  or  around  the 
suspended  soil  particles.  Ullrich's  model  can  be  used  in  a  computer  code 
when  the  soil  volumetric  strain  is  positive,  and  another  diffusion  model 
used  when  the  soil  volumetric  strain  is  negative. 
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4.  Morrison  (1979a) 

Morrison  reviewed  several  previous  studies  related  to  the  pore  air 
effect.  He  questioned  whether  Zernow's  vented  shock  tube  sample 
adequately  simulated  a  semi-infinite  soil  column,  but  did  not  attempt  to 
correct  Zernow’s  linear  solution  for  the  vented  case.  Instead,  he 
presented  a  derivation  of  the  nonlinear  equation  for  three-dimensional 
adiabatic  diffusion  of  air  through  a  rigid,  porous,  isotropic  medium, 
using  air  density  as  the  dependent  variable.  The  derivation  is  slightly 
simpler  when  pressure  is  used  as  the  dependent  variable,  and  is  presented 
in  Appendix  F,  including  a  pseudo-linear  form  simple  enough  for  even  hand 
calculation. 

Morrison  indicated  that  a  finite  difference  code  (presumably  one¬ 
dimensional)  had  been  developed,  in  which  the  porous  medium  behaves 
elastically  (and  presumably  linearly)  in  compression,  but  allows  air  flow 
between  zones  when  expanded.  Neither  equations  nor  numerical  results  were 
presented,  but  the  equations  of  motion  were  ascribed  to  the  WONDY  finite 
difference  code.  Soil  permeability  is  input  initially,  then  altered  as 
the  soil  porosity  changes.  The  computational  sequence  is  as  follows: 

a)  Using  previous  total  stress  and  total  density  as  inputs  to  WONDY, 
solve  equations  of  motion  to  obtain  new  displacements. 

b)  Calculate  new  strain  and  porosity. 

c)  Calculate  new  effective  stress. 

d)  Calculate  new  permeability. 

e)  Assuming  adiabatic  compression  without  flow,  calculate 
intermediate  pore  air  pressure. 

f)  Allowing  adiabatic  flow,  calculate  new  pore  air  pressure. 
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g)  Calculate  new  pore  air  density. 

h)  Calculate  new  total  stress. 

i)  Calculate  new  total  density. 

j)  Return  to  (a). 

5.  Morrison  (1979b) 

In  a  briefing  to  the  1979  DNA  Spall  Workshop,  Morrison  described  an 
experimental  device  designed  to  apply  transient  negative  (gage)  airblast 
pressure  to  the  top  surface  of  an  unjacketed,  vertical,  cylindrical  soil 
sample,  in  an  attempt  to  isolate  the  negative  phase  portion  of  the  pore 
air  effect.  Sixteen  millimeter  high  speed  photo  movies  were  shown 
illustrating  the  dramatic  influence  of  porosity  and  particle  gradation  on 
near-surface  particle  motion,  as  summarized  below: 


Soi  1 

Type 

Pore  Air  Flow 

Particle  Motion 

sand 

very  little 

particles  retain 
initial  relative 
position 

MISERS 

BLUFF 

[very  little] 

local  fluidization 
and  extensive  mixing 
of  layers 

gravel  with 
fines 

extensive 

fines  carried  through 
gravel  to  surface; 
little  motion  of  gravel 

After  reviewing  previous  data  on  airblast  penetration  from  [Hampton 
(1964)],  and  near-surface  soil  particle  lofting  from  [Zernow,  et  al 
(1973)],  Morrison  presented  additional  data  on  the  effect  of  sub- 
atmospheric  airblast  pressure  in  MISERS  BLUFF  1 1—2  to  show  pore  air 
expansion  effects.  A  definite  correlation  was  established  between 
negative  airblast  pressure  and  both  vertical  extensional  strain  and 
vertical  particle  velocity.  He  then  described  the  calculational  model  for 
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pore  air  expansion  discussed  in  [Morrison  (1979a)],  and  presented  some 
preliminary  results  in  which  peak  displacement  and  velocity  were  plotted 
as  functions  of  permeability  and  weapon  yield  for  a  given  peak  (negative) 
airblast  pressure.  Morrison  concluded  that  permeability  may  be  an 
important  site  characterization  factor. 

6.  Rosenblatt,  Qrphal  and  Hassig  (1979) 

In  another  presentation  to  the  1979  DNA  Spall  Workshop,  Rosenblatt 
et  al  reported  a  fundamental  approach  to  analysis  of  the  pore  air  effect, 
using  the  DICE  computer  code.  They  first  presented  a  concise  but 
comprehensive  summary  of  previous  observations  concerning  the  pore  air 
effect: 

a)  The  magnitude  of  near-surface  ground  response  to  the  small 
secondary  airolast  peak  (commonly  called  "repete" )  suggests  the  ground  has 
been  highly  dilated  and  has  a  relatively  low  impedance  (for  soil). 

b)  The  very  low  propagation  velocity  of  the  above  secondary 
compression  signal  is  consistent  with  the  assumption  of  a  two  phase  medium. 

c)  A  secondary  upward  velocity  appears  to  be  associated  with  the 
arrival  of  the  airblast  negative  phase,  and  is  hypothesized  to  be  caused 
by  upward  expansion  (flow)  of  pore  air. 

d)  Inclusion  of  Ullrich's  "no  flow"  adiabatic  expansion  model 
significantly  improves  the  correlation  between  calculated  and  measured 
near-surface  soil  vertical  particle  motion  during  the  negative  airblast 
phase,  for  dry,  high-porosity  soils. 

e)  Then  current  pore  air  expansion  models  made  no  provision  for  pore 
air  flow  or  the  associated  viscous  drag  or  seepage  force  on  the  soil 
skeleton.  Rosenblatt  et  al  suggested  that  the  pore  air  flow  mechanism  is 


APPLIED  RE/6RPCH  R//OCIRTE/,inC. 


C 


an  important  part  of  the  pore  air  effect. 

f)  The  shock  tube  experiments  by  [Zernow  et  al  (1973)]  showed 
conclusively  that,  under  their  test  conditions,  gas  permeation  of  the  sand 
column  was  necessary  to  produce  significant  upward  column  motion. 

g)  The  above  experiments  also  showed  that  upward  soil  particle 
velocities  are  enhanced  by:  higher  peak  surface  airblast  pressure,  more 
rapid  overpressure  decay,  smaller  particle  size,  and  an  underlying 
impermeable  layer;  and  are  reduced  by:  very  high  porosity  and 
permeability,  and  very  low  permeability.  [In  this  regard,  soil 
susceptibility  to  the  pore  air  effect  is  similar  to  soil  susceptibility  to 
frost  heave.] 

Several  questions  were  also  raised: 

a)  Is  the  pore  air  effect  two-dimensional,  i.e.,  are  horizontal 
motions  also  affected? 

b)  Is  the  negative  airblast  phase  required  to  produce  the  pore  air 
effect,  or  just  an  upward  pore  air  flow? 

c)  What  is  the  influence  of  soil  tensile  strength? 

d)  How  does  pore  air  related  behavior  scale  between  HE  and  NE 
surface  bursts? 

Specific  objectives  of  the  reported  effort,  which  was  just  beginning 
at  the  time  of  the  presentation,  were: 

a)  Formulate  a  one-dimensional  mathenatical  model  of  "forward"  and 
"reverse"  air  flow  into  the  ground,  and  "validate"  it  against  reverse 
percolation  laboratory  data  reported  by  [Zernow,  et  al  (1973)]. 

b)  Using  the  above  model,  calculate  relevant  vertical  ground  motions 
for  the  MISERS  BLUFF  II  experiments. 
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c)  Calculate  ground  motions  for  a  1  Iff  nuclear  surface  burst,  to 
assess  the  importance  of  two-phase  phenomena  for  nuclear  yields  of 
interest  and  to  investigate  scalability. 

Adaptation  of  the  DICE  code  to  handle  soil/air  interaction  was 
outlined,  but  the  definition  of  stress  in  the  soil  skeleton  appeared  to  be 
different  from  the  concept  of  effective  stress  used  in  conventional  soil 
mechanics.  The  formulation  did  allow  permeabi 1 ity  to  vary  with  soil 
porosity  over  a  wide  range  of  porosity,  as  well  as  with  soil  particle 
size.  The  results  of  this  work  have  since  been  presented  in  more  detail 
by  [Rosenblatt,  Hassig  and  Orphal  (1982)],  and  are  discussed  below. 

7.  Morrison,  Berglund  and  Kelly  (1979) 

Morrison,  Berglund  and  Kelly  performed  a  combined  experimental  and 
theoretical  study  of  soil  subjected  to  transient  negative  (gage)  airblast 
pressure.  Both  their  shock  tube  experiments  and  numerical  computer 
calculations  were  one-dimensional.  The  study  deliberately  provided  an 
isolated  view  of  that  portion  of  the  pore  air  effect  caused  by  the 
negative  phase  of  an  explosively-produced  air  shock.  The  main  problem 
with  such  an  approach  is  that  in  both  chemical  and  nuclear  explosions  a 
positive  airblast  phase  precedes  the  first  negative  phase,  so  the  soil  is 
not  apt  to  have  a  zero  internal  pore  air  pressure  distribution  when  the 
external  airblast  negative  pressure  phase  begins.  Thus,  the  positive  and 
negative  airblast  phase  effects  cannot  be  isolated  from  each  other.  Of 
course,  if  the  entire  process  were  linear,  the  positive  and  negative  phase 
effects  could  be  considered  separately  and  the  results  superposed. 

However,  there  is  good  reason  to  believe  that  not  all  aspects  of  the  pore 
air  phenomenon  are  linear,  particularly  the  soil  motions. 
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The  study  had  three  stated  primary  objectives: 

a)  form  a  basic  analytical  model  that  describes  soil  motion  due  to 
pore  air  expansion; 

b)  develop  an  experimental  apparatus  and  technique  that  can 
illustrate  and  measure  the  effects  of  pore  air  expansion;  and 

c)  provide  a  comparison  between  data  and  theory  so  theoretical 
limitations  and  possibilities  for  future  work  can  be  defined. 

Expansion  of  pore  air  due  to  pressure  differentials  with  depth,  and 
actual  pore  air  flow  throgh  the  soil  skeleton  are  treated  as  separate  soil 
lofting  mechanisms,  with  emphasis  on  the  former.  However,  an  effective 
stress  (or  soil  skeleton  oriented)  approach  shows  that  the  pore  air 
pressure  gradient  associated  with  flow  constitutes  a  distributed  body 
force  acting  throughout  the  soil  skeleton,  and  it  is  this  distributed 
seepage  force  which  helps  cause  lofting  or  spall.  Thus  pore  air  expansion 
and  pore  air  flow  are  not  two  phenomena,  but  slightly  different  aspects  of 
the  same  single  phenomenon. 

The  computer  code  employed  was  a  revision  of  the  one  described  by 
[Morrison  (1979a)],  and  employed  the  following  features: 

a)  a  one-dimensional  wave  equation,  using  total  stress  and  total 
density,  and  including  gravity  and  artificial  soil  viscosity; 

b)  a  soil  stress-strain  relation  using  effective  stress  but  no 
vi scosity; 

c)  Ullrich's  adiabatic  relation  between  pore  air  pressure,  soil 
porosity  and  soil  volumetric  strain.  Equation  (E.9); 

d)  a  total  mass  conservation  equation  which  assumes  no  pore  air 
flow;  and 
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e)  a  one-dimensional,  nonlinear  diffusion  equation  for  isothermal 
pore  air  flow  through  a  rigid,  porous  soil  skeleton.  Equation  (A. 5),  which 
assumes  the  validity  of  Darcy's  law.  Equation  (A.l). 

Three  different  expressions  are  given  for  the  viscous  stress  used  in 
the  one-dimensional  wave  equation:  one  in  the  text,  one  in  an  appendix, 
and  one  in  the  computer  program  listing.  No  mention  is  made  of  the  pore 
air  seepage  force,  or  of  the  fact  that  both  adiabatic  and  isothermal 
conditions  are  assumed  in  the  same  set  of  equations.  Therefore, 

Appendix  F  of  this  report  presents  the  derivation  of  a  pseudo-linear 
diffusion  equation  for  adiabatic  pore  air  flow  through  a  rigid,  porous 
soil  skeleton  which  is  even  simple  enough  for  hand  calculation. 

The  only  mention  of  bottom  boundary  conditions  imposed  in  the  above 
calculations  is  the  observation  that  "Dramatic  differences  [50  to  75  mm 
(2  to  3  in)]  in  deflections  were  observed  when  the  total  sample  depth  was 
varied  by  only  25  mm  (1  in)."  However,  only  two  tests  are  reported  in 
which  only  the  sample  depth  was  varied.  Mention  is  made  of  a  single 
calculation  having  been  made  using  an  airblast  input  with  both  positive 
and  negative  phases,  but  the  results  are  not  reported.  A  theoretical 
analysis  of  the  linear  pore  pressure  response  in  a  rigid,  porous  solid, 
subjected  to  surface  airblast  loading  having  both  a  positive  and  a 
negative  phase  is  presented  in  Appendix  G  of  this  report. 

8.  Merkle  (1980) 

Merkle  examined  current  soil  spall  theories,  including  Ullrich's  pore 
air  expansion  mechanism.  The  pore  air  expansion  equations  were  derived 
using  conventional  soil  mechanics  nomenclature,  and  were  shown  not  to  be 
restricted  to  one-dimensional  effects.  They  define  the  propagation 
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velocity  of  a  pressure  wave  in  a  two-phase  fluid,  in  which  one  phase  is 
compressible  but  massless,  while  the  other  phase  has  mass  but  is 
incompressible. 

9.  Labreche  (1980) 

Labreche  described  a  series  of  finite  difference  computer 
calculations,  using  the  codes  WONDY  and  WONDY/POREAIR,  to  study  the 
differences,  if  any,  between  pore  air  effects  for  a  1  MT  nuclear  loading 
and  a  high  explosive  loading  scaled  to  1  MT,  both  with  the  same  peak 
side-on  airblast  overpressure.  The  HE  loadings  selected  for  comparison 
were  from  MISERS  BLUFF  Phase  I  Event  2  and  Phase  II  Event  1.  The  soil 
profiles  used  above  were  adjusted  and  scaled  versions  of  the  MISERS  BLUFF 
Phase  I  and  Phase  II  test  site  profiles.  Several  other  calculations  were 
made  in  support  of  vertical  shock  tube  studies,  and  for  comparison  with 
results  obtained  by  others  using  a  no  flow  model.  The  soil  total  stress- 
strain  model  used  in  the  WONDY  code  was  the  AFWL  Engineering  Model. 

Two  sets  of  calculations  were  run  for  the  HE/NE  comparison  cases: 
one  set  used  the  code  WONDY,  and  ignored  the  pore  air  effect;  the  other 
used  the  code  WONDY/POREAIR  and  included  the  pore  air  effect.  The 
difference  between  the  soil  particle  motions  predicted  by  WONDY/POREAIR 
and  by  WONDY  alone  was  taken  as  the  measure  of  the  pore  air  effect.  When 
the  WONDY/POREAIR  code  was  used  the  pore  air  diffusion  equations  in 
POREAIR  were  bypassed  when  the  soil  extensional  strain  reached  a 
predetermined  limiting  value,  determined  by  the  maximum  staDle  soil 
porosity.  At  this  point  the  soil  particles  were  assumed  to  be  no  longer 
in  contact,  and  to  undergo  ballistic  motion.  At  this  point  the 
airblast-loaded  surface  was  effectively  shifted  down  to  the  top  of  the 
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first  nonspalled  element.  During  initial  compression,  WONDY/POREAIR 
assumed  soil  porosity,  permeability,  pore  air  density  and  pore  air 
pressure  to  all  be  constant.  Numerical  stability  was  studied  by  varying 
the  WONDY  propagation  time  step,  at,  and  the  necessarily  smaller  POREAIR 
diffusion  time  step,  a t/N. 

10.  Morrison,  et  al  (1981) 

Morrison,  Labreche,  and  Lamb  reported  in  detail  the  results  of 
evacuation  chamber  tests  and  HE/NE  comparison  calculations  which  had  been 
outlined  by  [Labreche  (1980)].  They  also  gave  preliminary  results  of 
vertical  shock  tube  tests,  in  which  the  airblast  load  had  both  a  positive 
and  a  negative  phase.  They  concluded  that  pore  air  effects  are  sensitive 
to  soil  permeability,  which  in  turn  is  affected  by  soil  particle  gradation 
and  degree  of  saturation,  as  well  as  to  the  water  phase  in  a  partly 
saturated  soil.  Saturated  soil  exhibits  no  pore  air  effects,  because 
there  is  no  connected  gaseous  pore  air  phase. 

In  this  report,  pore  air  expansion  of  soil  is  defined  as  "expansion 
of  the  soil  which  results  when  the  pore  air  pressure  in  the  soil  voids 
exceeds  the  confining  stress  of  the  soil  matrix".  A  literal 
interpretation  of  the  above  definition  implies  a  negative  effective 
stress.  This  was  not  intended.  The  situation  envisioned  by  the 
definition  is  that  of  a  finite  soil  mass,  in  which  the  pore  air  pressure 
within  some  interior  region  exceeds  the  pore  air  pressure  over  some 
surface  region,  thus  producing  an  outward  pore  air  flow,  the  seepage  force 
from  which  tends  to  cause  soil  particle  contact  loss. 

Near-surface  displacement-time  curves  from  evacuation  chamber  tests 
were  obtained  from  high  speed  photographs  and  also  double  integration  of 
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accelerometer  data.  The  Integrated  accelerometer  records  indicated  much 
larger  peak  upward  displacements  than  did  the  photographs,  and  the 
accelerometers  came  to  rest  nearer  the  surface  than  they  had  been  placed. 
This  led  to  the  conclusion  that  the  near-surface  accelerometers  did  not 
move  with  the  surrounding  soil  (perhaps  at  least  in  part  because  their  air 
drag  characteristics  were  different  from  those  of  soil  particles). 

11.  Rosenblatt,  Hassig  and  Orphal  (1982) 

Rosenblatt,  Hassig  and  Orphal  described  some  of  the  theoretical 
fundamentals  used  to  modify  the  DICE  code  to  treat  pore  air  effects 
related  to  lofting,  and  gave  results  of  their  calculations.  This  is  the 
work  outlined  in  their  1979  DNA  Spall  Workshop  presentation  discussed 
above. 

The  conventional  soil  mechanics  definition  of  effective  stress  is 
presented  in  the  theoretical  discussion,  but  some  of  the  calculations 
(e.g.  Figure  4.10(a))  appear  to  show  a  large  effective  stress  at  the 
unjacketed  surface  of  a  soil  mass  loaded  only  by  pore  air  pressure.  Since 
the  effective  stress  at  such  a  surface  is  always  zero,  it  is  not  clear 
that  the  DICE  code  actually  used  the  conventional  soil  mechanics 
definition  of  effective  stress.  This  situation  is  encountered  in 
conventional  geotechnical  engineering  practice  in  connection  with  a  change 
of  water  level  in  a  reservoir  or  river.  It  also  happens  continually  in 
tidal  areas. 

Results  of  the  DICE  code  calculations  showed  reasonable  agreement 
with  some  of  Zernow's  shock  tube  test  results,  although  that  was  not  the 
stated  objective  of  the  calculations.  In  fact,  the  calculations  used 
different  bottom  displacement  and  pore  pressure  boundary  conditions  than 
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those  appropriate  for  Zernow’s  samples;  they  also  kept  the  airblast  decay 
rate  constant  and  varied  soil  permeability  in  obtaining  results  which  were 
compared  with  those  of  tests  in  which  the  airblast  decay  rate  varied  and 
the  soil  permeability  remained  constant.  Nevertheless,  the  DICE  code 
results'  qualitative  agreement  with  Zernow' s  indicate  that  the  code  has 
the  capability  to  model  the  essential  features  of  the  pore  air  effect. 

C.  Simple  Discrete  Models 

The  above  pore  air  effect  analyses  employ  sets  of  equations  which  are 
only  partially  coupled,  in  the  sense  that  the  equations  describing  stress 
wave  propagation  and  those  describing  pore  air  diffusion  are  solved 
sequentially,  rather  than  simultaneously  for  each  time  increment.  The 
difficulty  of  formulating  and  solving  a  completely  coupled  set  of 
equations  for  the  pore  air  effect  was  recognized  early  in  this  study. 
Therefore,  to  gain  insight  into  both  the  physical  mechanisms  and  the 
mathematical  processes  involved,  some  simple  mathematical  models  were 
constructed,  the  behavior  of  which  could  be  analyzed  in  closed  form. 

The  first  model,  analyzed  in  Appendix  H,  actually  arose  during  a 
discussion  of  the  response  of  an  unjacketed  soil  test  sample  immersed  in 
water  in  a  closed  pressure  vessel,  when  the  fluid  pressure  suddenly 
increased.  This  of  course  is  the  reservoir  or  tidal  basin  problem 
mentioned  previously.  The  simple  piston  model  analyzed  in  Appendix  H  is 
massless,  so  there  are  no  inertial  effects.  The  main  thing  this  model 
shows  is  that  when  a  saturated,  un jacketed  soil  sample  is  subjected  to  a 
sudden  increase  in  external  fluid  pressure,  the  sample  first  undergoes 
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undrained  composite  compression,  in  which  the  soil  skeleton  and  the 
internal  pore  fluid  experience  the  same  volume  decrease.  Subsequent 
drainage  allows  compressed  pore  fluid  to  flow  into  the  soil  skeleton,  and 
the  soil  skeleton  to  expand  to  its  original  volume  (if  elastic). 

The  second  model,  analyzed  in  Appendix  I,  is  the  same  as  the  first 
except  that  both  the  soil  skeleton  and  the  pore  air  have  mass.  With  this 
model,  initial,  instantaneous  undrained  deformation  does  not  occur  (see 
Equation  (1.44)),  because  of  the  inertia  of  both  soil  skeleton  and  pore 
fluid.  However,  the  numerical  evaluation  of  Equation  (1.62)  shows  that  a 
small  deformation  of  the  soil  skeleton,  associated  with  pore  air 
compression,  does  occur  very  rapidly.  Subsequent  deformation  of  the  soil 
skeleton,  associated  with  pore  air  flow,  occurs  much  more  slowly.  The 
final  equilibrium  condition  of  this  model  is  the  same  as  that  for  the 
massless  model.  The  internal  pore  fluid  eventually  attains  the  same 
pressure  as  the  external  fluid,  and  the  soil  skeleton  returns  to  its 
initial  configuration. 

Since  both  the  above  models  are  linear,  their  response  to  any 
prescribed  airblast  input  can  be  obtained  by  superposition. 

D.  General  Equations 

The  advantage  of  the  above  two  discrete  models  of  soil  pore  air 
behavior  is  that  they  are  relatively  simple  and  linear,  so  that  their 
response  can  be  obtained  in  closed  form.  This  permits  parametric  studies 
and  physical  interpretation  of  the  results.  Nevertheless  what  is  obtained 
is  an  exact  solution  to  an  approximate  problem,  rather  than  an  approximate 
solution  to  the  real  problem.  The  general  equations  describing  the  "real" 
one-dimensional  pore  air  problem  are  developed  in  Appendix  M.  They  are  a 
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set  of  coupled,  nonlinear,  second  order  partial  differential  equations 
involving  soil  skeleton  displacement  and  pore  air  pressure.  The  most 
practical  approach  to  their  solution  appears  to  be  by  the  method  of  finite 
differences. 
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SECTION  V 


Prediction  of  Negative  Airblast  Overpressures 
From  Near-Surface  Explosions 

It  will  be  shown  in  Section  VI  that  a  major  portion  of  the  shallow 
spall  observed  in  experiments  with  near-surface  explosions,  detonated 
over  granular  materials  with  essentially  zero  tensile  strength,  results 
from  the  negative  overpressure  portion  of  the  airblast  loading.  Thus, 
the  negative  airblast  overpressure  must  be  known  in  order  to  develop  a 
model  to  predict  this  phenomenon.  A  review  of  the  literature  revealed 
that  no  method  was  available  to  accurately  define  the  negative  airblast 
overpressure  from  near-surface  explosions.  Accordingly,  available  data 
and  calculations  were  reviewed,  and  a  procedure  developed  for  predicting 
the  maximum  negative  airblast  overpressure  as  a  function  of  scaled  range. 

Figure  5.1  presents  measured  maximum  negative  airblast  overpressure, 
aP~,  versus  scaled  range,  R-,  from  several  high  explosive  near-surface 
detonations.  Note  that  the  maximum  measured  value  of  aP~  is 
approximately  48  kPa.  Lack  of  data  at  higher  values  of  aP”  indicates 
the  difficulty  of  developing  instrumentation  to  make  such  measurements. 

Figure  5.1  also  presents  the  results  of  numerical  calculations 
[Needham,  1969)]  for  a  500  ton  high  explosive  surface  tangent  event.  The 
calculated  values  for  30  kPa  and  below  are  in  reasonable  agreement  with 
the  measured  data  in  Figure  5.1,  but  indicate  a  somewhat  flatter  decay 
rate  with  scaled  range  than  an  upper  bound  to  the  measured  data.  The 
calculations  for  surface  tangent  and  slightly  aboveground  detonations 
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indicate  the  formation  of  a  toe  of  detonation  products  near  ground  level, 
which  signif icantly  influences  the  airblast  parameters.  The  value  of 
"R  a  0.10  km/(kt)l/3,  at  which  a  jump  in  aP-  occurs,  is  the  scaled 
range  at  which  a  normal  shock  is  formed.  At  smaller  values  of  R",  the 
results  approach  a  maximum  value  of  -1  atmosphere. 

An  analytical  approximation  to  both  the  measured  and  calculated  data, 
consisting  of  two  straight  lines,  is  also  shown  in  Figure  5.1.  The 
equations  for  this  proposed  approximation  are: 

aP-  =  Se.OCR]"0-1^,  |<Pa  for  R<0.15  km/(kt)l/3  (5.1) 

and 

aP~  =  3.98[R]-1*329f  |<Pa  for  R>0.15  km/(kt )1/3  (5.2) 

This  analytical  approximation  is  nearly  an  upper  bound  to  the  measured 
points.  The  breakover  point  and  slope  of  the  upper  portion  of  the 
analytical  approximation  are  based  on  engineering  judgment,  and  are  not 
in  exact  agreement  with  the  calculated  values. 

Measured  airblast  results  from  several  large  yield  (lOkt  -  lOmt) 
nuclear  events  conducted  in  the  Pacific  are  summarized  in  scaled  format 
in  Table  5.1.  The  aP“  versus  "R  nuclear  data  are  plotted  in 

Figure  5.2.  Again  there  are  no  data  at  the  larger  aP-  values,  i.e., 

for  aP~>28  kPa.  Results  of  calculations  utilizing  the  1  kt  nuclear 
standard  [Needham  (1975)]  are  also  presented  in  Figure  5.2  for 
comparison.  The  calculated  values  again  indicate  a  flatter  decay  rate 
than  an  upper  bound  to  the  measured  data,  but  in  this  instance  the 
calculated  values  are  larger  at  most  scaled  ranges.  The  calculations 
indicate  a  breakover  point  at  an  K  of  approximately  0.12  km/(kt)l/3. 

At  smaller  values  of  "R  the  results  approach  a  maximum  value  of 
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approximately  -1/3  atmosphere  instead  of  the  -1  atmosphere  indicated  for 
high  explosives. 

The  high  explosive  approximation  for  the  peak  negative  airblast 
overpressure  was  corrected  by  a  factor  of  2W  to  account  for  radiation 
loss  occurring  in  a  nuclear  detonation  (1  kt  nuclear  =  1/2  kt  high 
explosive).  This  corrected  line  is  shown  in  Figure  5.2.  It  provides  an 
excellent  fit  as  an  upper  bound  to  the  measured  nuclear  data.  Thus,  it 
appears  that  a  single  straight  line,  corrected  by  2W,  can  be  used  to  fit 
both  the  nuclear  and  high  explosive  measured  data  at  larger  scaled 
ranges.  The  nuclear  breakover  point  is  maintained  at  0.15  km/(kt)^, 
and  the  slope  of  the  proposed  upper  line  in  the  nuclear  approximation  is 
reduced  somewhat  from  the  high  explosive  approximation.  The  proposed 
equations  for  the  nuclear  approximations  are: 

aP-  =  27.3[R]-°-132,  kPa  for  R<0.15  km/(kt)1/3  (5.3) 

and 

aP-  =  2.90[R]-1-313,  kPa  for  R>0.15  km/(kt)l/3  (5.4) 

The  two  analytic  approximations  for  predicting  peak  negative  airblast 
overpressure  versus  scaled  range  are  presented  in  Figure  5.3.  The  yield 
for  the  appropriate  type  of  explosive  is  utilized  without  correction  for 
this  figure.  The  curves  provide  upper  bounds  to  the  available  measured 
data.  Considerable  uncertainty  exists  at  scaled  ranges  less  than 
0.15  km/(kt)^,  and  the  above  results  should  be  used  with  caution  for 
scaled  ranges  below  that  value. 
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SECTION  VI 


Review  of  Spall  Data  From  HE  and  NE  Detonations 

Spall  has  been  recently  studied  by  [Stump  and  Reinke  (1980),  Merkle 
(1980),  and  Auld,  et  al  (1981)].  For  the  purposes  of  this  report,  the 
following  criteria  will  be  utilized  to  identify  spall  in  ground  motion 
records: 

1.  -lg  (-0.5  to  -2.0)  vertical  acceleration  dwell  (j>5  ms),  which  is 
identifiable  directly  on  a  vertical  acceleration  record,  or  as 
the  slope  of  a  vertical  velocity  record; 

2.  identifiable  impulsive  rejoin  signal  on  both  the  vertical  and 
horizontal  acceleration  records; 

3.  rejoin  amplitude  (>0.05  m/s)  observed  on  the  vertical  velocity 
record. 

These  criteria  are  similar  to  those  utilized  by  previous  investigators, 
but  with  numerical  values  specified. 

[Auld,  et  al  (1981)]  suggested  that  the  zone  of  spalled  material 
surrounding  a  near-surface  explosion  can  be  considered  to  consist  of  two 
parts:  a  bowl-shaped  volume  of  material,  designated  the  "coupled  spall 
region";  and  a  shallow  extension  of  the  basic  spalled  volume,  designated 
the  "negative  airblast  wing  region".  Figure  6.1  illustrates  these  two 
zones  as  they  might  occur  from  a  typical  near-surface  detonation.  Any  of 
the  possible  spall  mechanisms,  or  a  combination  of  these  mechanisms  may 
be  responsible  for  the  spall  observed  in  the  coupled  region,  e.g.,  direct 
waves,  head  waves,  reflected  waves,  or  surface  waves.  The  predominant 
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spall  mechanism  in  the  negative  airblast  wing  region  is  the  negative 
overpressure  portion  of  the  airblast  loading.  Spall  can  occur  in  this 
region  only  when  the  applied  peak  negative  airblast  overpressure  exceeds 
the  tensile  strength  of  the  near-surface  material.  Accordingly,  many 
tests  do  not  exhibit  the  negative  airblast  spall  wing  because  the 
explosive  charge  was  slightly  buried  (thereby  suppressing  the  airblast), 
the  surface  material  had  a  significant  tensile  strength  in  comparison  to 
the  applied  peak  negative  airblast  overpressure,  or  the  instrumentation 
was  placed  in  holes  backfilled  with  grout  instead  of  with  negligible 
tensile  strength  material,  e.g.,  sand.  Note  that  "soil  matching"  grout 
can  be  utilized  to  backfill  instrumentation  holes  in  the  coupled  spall 
zone  and  spall  can  still  be  detected. 

[Auld,  et  al  (1981)]  also  suggested  that  a  simple  model  can  be  used 
to  estimate  the  depth  of  spall  in  the  negative  airblast  wing  for  a 
material  with  zero  tensile  strength.  This  model  equates  the  geostatic 
vertical  total  stress  at  the  depth  of  spall  to  the  value  of  the  peak 
negative  airblast  overpressure  at  the  range  of  interest,  thereby 
obtaining  the  equation: 


zs  =  tP-h  (6.1) 

where 

z$  =  depth  of  spall  at  the  range  of  interest  (see  Figure  6.1); 

&P"  =  peak  negative  airblast  overpressure; 
y  =  soil  total  unit  weight. 

3  _ 

For  p  s  1900  kg/m  ,  and  aP  in  kPa,  Equation  6.1  reduces  to: 

Zs  =  (9.80665^(1900)  =  0,064 
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The  pea!'  negative  airblast  overpressure  must  be  known  to  calculate 
from  Equation  6.2.  Accordingly,  the  technique  for  predicting  peak 
negative  airblast  overpressure  developed  in  Section  V  should  be  utilized. 

Table  6.1  summarizes  the  high  explosive  events  which  were  examined  in 
detail,  to  develop  an  understanding  of  spall  phenomena  and  to  develop 
techniques  for  predicting  the  extent  of  spall  associated  with  a 
near-surface  detonation.  Emphasis  was  placed  on  large  yield  events 
(W>100  tons)  and  on  events  with  sufficient  instrumentation  to  produce  a 
good  definition  of  the  spalled  region.  The  events  summarized  in 
Table  6.1  were  conducted  on  a  wide  range  of  geologic  profiles,  ranging 
from  deep  dry  soil  to  layered  soil  over  rock,  with  and  without  high  water 
tables.  In  addition,  the  near-surface  materials  ranged  from  dry  powdery 
silt,  with  essentially  zero  tensile  strength,  to  highly  competent 
granite.  High  explosive  yields  ranged  from  a  fraction  of  a  ton  to  500 
tons.  Charge  configurations  consisted  of  slight  height  of  burst  (HOB), 
surface  tangent  spheres  (STS),  surface  tangent  cylinders  with 
hemispherical  caps  (STC),  half  buried  spheres  (HBS),  berms,  and  slight 
depth  of  burial  (DOB).  All  high  explosive  events  studied  had 
identifiable  regions  of  spall,  except  the  three  shots  conducted  on 
granite,  i.e.,  MINERAL  ROCK,  MINE  ORE,  and  MINE  LINDER.  There  may  have 
been  spall  associated  with  the  granite  events;  however,  it  is  reasonable 
to  assume  that  high  rock  strength  limited  the  spall  region  to  the  extent 
that  violent  cratering  motions  masked  its  presence. 

Spall  was  detected  on  only  one  nuclear  event,  JANGLE-U.  Ground 
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motion  records  were  also  examined  for  MIKE,  CACTUS,  KOA,  PRISCILLA, 
JOHNNIE  BOY,  SMALL  BOY,  and  JANGLE-S,  without  detecting  motions 
satisfying  the  previously  specified  spall  criteria.  The  high  water  table 
associated  with  events  conducted  in  the  Pacific  may  have  precluded 
spall.  Grouted  instrumentation  holes,  old  style  instrumentation,  and 
less  sophisticated  data  processing  techniques  may  have  contributed  to 
masking  of  spall  in  nuclear  events  conducted  at  the  Nevada  Test  Site 
(NTS).  In  addition,  the  ground  motion  instrumentation  was  placed  too 
deep  (e.g.,  1.5  m)  to  record  much  of  the  spall  phenomena  on  all  of  the 
nuclear  events.  Figure  7.1  shows  the  extent  to  which  spall  in  the 
negative  airblast  wing  region  can  be  overlooked  by  using  a  minimum  gauge 
depth  of  this  magnitude.  Figure  7.1  also  indicates  that  using  a  much 
shallower  minimum  gauge  depth  for  high  explosive  experiments  increases 
the  probability  of  detecting  spall.  Note  that  high  explosive  events  at 
NTS  have  produced  significant  spall  regions.  However,  it  should  not  be 
concluded  on  the  basis  of  these  results  that  nuclear  events  do  not 
produce  spall. 

The  extent  of  the  spalled  volume  of  material  can  be  defined  by  two 
quantities:  radius  of  spall,  R  ,  and  maximum  depth  of  spall,  z$. 

Rs  will  vary  widely,  depending  on  whether  the  negative  airblast  spall 
wing  is  present  (see  Figure  6.1).  The  depth  of  the  negative  airblast 
spall  wing  is  a  slowly  decaying  function  of  range,  and  the  radius  of 
spall  for  zero  depth  can  only  be  crudely  approximated.  Accordingly,  R$ 
is  arbitrarily  defined  for  a  depth  of  0.5  m  in  this  report,  and  values 
corresponding  to  this  definition  are  given  in  Table  6.1.  The  radius  of 
spall  for  any  other  chosen  depth  can  be  estimated  by  using  Equation  6.2, 
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or  other  techniques  given  below.  In  the  coupled  spall  region  there  is 
essentially  no  difference  in  the  radius  of  spall  for  depths  of  zero  and 
0.5  m.  There  have  been  a  very  limited  number  of  close-in  measurements 
made  on  near- surface  detonations,  and  the  estimates  of  z$  are  based 
almost  entirely  on  extrapolation  of  data  usi.,g  a  preconceived  concept  of 
the  general  shape  of  the  spalled  volume.  Therefore,  there  is  a  great 
deal  of  uncertainty  associated  with  zg. 

Two  general  methods  were  utilized  to  estimate  the  values  of  R$  and 
zs  presented  in  Table  6.1.  The  first  method  consisted  of  evaluating 
ground  motion  records  using  the  established  spall  criteria.  Each 
instrumentation  location  was  then  labelled  as  showing  spall  (S), 
questionable  spall  (?S),  or  no  spall  (NS).  If  the  instrumentation  array 
is  sufficiently  large  and  dense,  the  zone  of  spalled  material  can  then  be 
estimated  as  shown  in  Figure  6.2.  A  second  method  was  developed  to 
estimate  the  radius  of  spall  from  the  magnitude  of  rejoin  amplitudes 
observed  on  vertical  velocity  records.  It  was  noted  that  the  radius  at 
which  the  near-surface  vertical  velocity  rejoin  amplitude  approached  zero 
correlated  well  with  the  radius  of  spall  defined  by  the  first  method. 

For  example,  vertical  velocity  rejoin  data  from  two  of  the  PRE-MINE 
THROW-I V  (PMT-IV)  events  are  presented  in  Figure  6.3,  for  a  depth  of 
0.23  m.  The  radius  of  spall  can  easily  be  estimated  by  extrapolating  the 
data  to  the  range  where  the  rejoin  amplitude  is  zero.  However,  this 
radius  of  spall  does  not  meet  the  above  definition  of  R$  (i.e., 
z  =  0.50  m);  and  since  the  PMT  events  had  airblast  spall  wings,  this 
slight  depth  difference  may  produce  significant  variations  in  the  radius 
of  spall.  Also  note  that  there  is  a  dip  in  the  rejoin  amplitude  data, 
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which  extrapolates  to  a  range  of  approximately  10  m.  This  range 
corresponds  well  to  the  R$  obtained  from  experiments  where  there  was  no 
negative  airblast  spall  wing,  and  suggests  that  high  quality  rejoin 
amplitude  versus  range  data  can  be  used  to  estimate  the  radius  of  spall 
for  both  coupled  and  negative  airblast  wing  regions. 

Figure  6.4  shows  rejoin  amplitude  plotted  against  scaled  range, 
R(m/0.5  ton  1  ),  for  all  PMT  events.  This  figure  illustrates  the  main 
difficulty  in  comparing  spall  data  from  different  events,  when  the  ground 
motion  instrumentation  for  each  event  is  placed  at  a  single  depth  which 
varies  from  event  to  event. 

These  data  suggest  that  the  scaled  value  of  FTs  for  z  =  0.50  m 
adequately  represent  all  PMT  events,  and  that  apparent  differences  in 
result  from  the  depth  of  the  ground  motion  instrumentation  rather 
than  from  differences  in  spall  phenomena.  The  insert  in  Figure  6.4 
illustrates  a  technique  for  crudely  estimating  the  in-situ  tensile 
strength  of  the  near-surface  material.  Equation  6.2  was  used  in 
conjunction  with  measured  values  of  &P~  to  estimate  the  maximum  depth 
of  spall  at  various  ranges,  for  a  zero  tensile  strength  material.  The 
difference  between  the  geostatic  vertical  total  stress  at  spall  depth  and 
the  applied  value  of  aP~  was  assumed  to  result  from  the  tensile 
strength  of  the  material  at  that  depth.  The  four  points  obtained  produce 
a  oj  versus  depth  curve  that  is  quite  plausible  for  that  particular 
site. 

The  values  of  spall  radius,  R$,  at  a  depth  of  0.5  m  obtained  by 
either  of  the  general  methods  previously  described  are  summarized  in 
Figure  6.5  by  plotting  scaled  radius  of  spall,  ^(m/lOO  ton^/^). 
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versus  depth  to  the  second  layer,  which  may  be  either  the  water  table  or 
competent  rock.  The  data  generally  fall  into  two  groups:  those  with  a 
negative  airblast  spall  wing,  and  those  without.  Yield  cube  root 
(W1'3)  scaling  for  the  radius  of  spall  appears  to  order  the  data. 

Neither  charge  configuration  nor  details  of  the  deeper  geologic  profile 
greatly  influences  Rs  when  the  negative  airblast  spall  wing  is 
present.  Of  course  the  near-surface  material  must  be  very  weak  to 
produce  a  negative  airblast  spall  wing.  R"s  is  influenced  by  the  charge 
configuration  when  the  negative  airblast  spall  wing  is  not  present,  and 
the  value  of  R$  increases  with  increased  charge  coupling 
(HOB-STS-STC-HBS-Berm) .  A  surface  tangent  cylinder  appears  to  be 
essentially  equivalent  to  a  half  buried  sphere  instead  of  a  surface 
tangent  sphere.  The  widest  data  scatter  occurs  for  the  STS  events. 
Similar  observations  have  previously  been  noted  with  regard  to  cratering 
parameters,  such  as  the  crater  radius  or  depth.  Again,  there  is  no 
indication  of  dependence  on  geological  features,  such  as  the  depth  to  the 
second  layer. 

R^s  for  the  JANGLE-U  event  is  also  shown  in  Figure  6.5.  This  event 
had  a  slight  depth  of  burial  (5.18  m)  and  should  have  produced  results 
comparable  to  the  bermed  high  explosive  event.  A  comparison  of  these  two 
data  points  will  be  used  in  Section  VII  to  estimate  a  nuclear  to  high 
explosive  spall  efficiency  factor,  and  thereby  develop  a  nuclear  spall 
prediction  technique. 

The  range  of  R”  values  for  the  various  charge  configurations  is 
summarized  in  Table  6.2  for  both  high  explosive  and  nuclear  events. 

The  dependence  of  Rs  on  charge  configuration  for  the  events  with  no 
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negative  airblast  wing,  i.e.,  with  only  a  coupled  spall  region,  suggests 
that  there  might  be  a  correlation  between  unsealed  radius,  Rg,  and 
crater  size.  However,  yield  cube  root  scaling  collapsed  the  data  better 
than  scaling  by  either  crater  radius  or  crater  volume. 

The  relation  between  maximum  depth  of  spall,  zs,  and  depth  to  the 
second  layer  (either  the  water  table  or  competent  rock)  is  shown  in 
Figure  6.6.  It  appears  that  there  is  no  strong  dependence  on  charge 
conf iguration,  as  was  the  case  in  Figure  6.5.  However,  a  dependence  on 
geology  is  indicated.  The  maximum  depth  of  spall  appears  to  be  limited 
by  the  depth  to  the  water  table  when  there  is  no  shallower  rock  layer. 
This  limitation  is  indicated  by  the  data  for  which  the  spall  zone  ends  in 
soil  falling  on  or  to  the  left  of  the  line  for  which  maximum  depth  of 
spall  equals  depth  to  the  second  layer.  However,  as  indicated  by  the  two 
data  points  for  which  the  spall  zone  ends  in  rock,  when  spall  does 
penetrate  into  weak  to  intermediate  strength  rock  (sandstone  or  shale), 
the  water  table  may  not  stop  it.  The  data  are  limited  but  suggest  that 
the  maximum  depth  of  spall  is  less  for  the  case  of  a  water  table  over  a 
rock  layer  than  for  the  case  of  dry  soil  extending  to  the  rock  layer. 

The  relation  between  maximum  depth  of  spall  and  yield  is  shown  in 
Figure  6.7.  A  clear  dependence  on  yield  is  indicated.  In  addition,  the 
geology  dependence  previously  seen  in  Figure  6.6  is  again  evident,  and 
there  is  a  suggestion  that  charge  configuration  may  also  be  of  some,  but 
lesser  importance.  These  data  are  insufficient  to  account  for  detailed 
differences  in  geology,  but  the  following  equations  can  be  utilized  to 
estimate  z$  for  a  surface  tangent  charge  configuration  (either  STC  or 
STS): 
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(# 

z$  =  depth  of  water  table* 
or  3.2  (W,tons)*^,m 

®  whichever  is  smaller.  (6.3) 

A  similar  equation  applies  to  the  case  where  the  charge  is  half  buried 
(HBS) : 

®  2$  =  depth  of  water  table* 

or  5.0(W,tons)^,m.  (6.4) 

Neither  Equation  (6.3)  nor  Equation  (6.4)  applies  when  a  competent  rock 
^  layer  is  encountered  above  2s. 


*If  there  is  a  rock  layer  at  a  depth  of  2S  or  less  disregard  the  water 
table. 
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SECTION  VII 


Prediction  of  Spall  Region  for  Near-Surface  Detonations 
A.  Introduction 

The  prediction  techniques  presented  in  this  section  are  based  on  the 
negative  airblast  predictions  from  Section  V  and  the  measured  spall  data 
presented  in  Section  VI.  Note  that  spall  was  detected  on  only  one 
nuclear  event.  Therefore,  the  extrapolation  from  high  explosive  to 
nuclear  involves  a  great  deal  of  uncertainty,  and  these  results  should  be 
used  with  caution.  Also,  no  spall  was  observed  in  very  competent  rock 
such  as  granite. 

As  previously  discussed,  the  volume  of  spalled  material  from  a 
near-surface  detonation  can  be  characterized  by  two  parameters,  the 
maximum  radius  of  spall,  Rs,  and  the  maximum  depth  of  spall,  z$. 

Rs  is  arbitrarily  defined  to  be  at  a  depth  of  0.5  m,  and  may  be 
associated  with  either  the  coupled  spall  region  or  the  negative  airblast 
wing,  when  this  region  is  present  {see  Figure  6.1). 

To  make  predictions  of  R$  and  z$,  the  following  information  must 
be  known: 

1.  Yield  (W):  tons  of  TNT  or  kt  of  nuclear; 

2.  Charge  configuration:  surface  tangent  or  half  buried  sphere 

(small  HOB  for  nuclear  is  assumed); 

3.  Geology 

a.  depth  to  water  table 

b.  depth  to  rock  layer  and  type  of  rock 

c.  nature  of  the  near-surface  material  (does  it  have  a 
negligible  tensile  strength?,  i.e.,  oy<10  kPa?^ 


51 


APPLIED  RE/6APCH  A//OCIATE/,mC. 


B.  Spall  Associated  With  High  Explosive  Detonations 

When  the  near-surface  material  has  a  negligible  tensile  strength  and 
airblast  is  present,  a  negative  airblast  wing  region  can  be  expected.  In 
this  instance  the  expected  scaled  maximum  radius  of  spall  can  be 
estimated  from  Figure  6.5  to  be 

IT  = - =  32  m/ton1/3  (7.1) 

s  (ioor/J 

Additional  definition  of  the  shape  of  the  negative  airblast  wing  region 

can  be  obtained  through  the  use  of  Figure  7.1.  This  figure  was  derived 

by  combining  the  predictions  of  aP-  from  Figure  5.3  with  Equation  6.2. 

Since  the  tensile  strength  of  the  near-surface  material  is  assumed  to  be 

zero,  the  results  will  be  conservative.  The  predicted  depth  of  spall  can 

be  easily  obtained  for  any  scaled  radius  from  Figure  7.1,  and  the  shape 

of  the  negative  airblast  wing  region  estimated.  Note  that  the 

theoretical  scaled  radius  of  HE  spall  for  a  depth  of  0.5  m  shown  in 

Figure  7.1  is  0.520  km/kt^3.  This  scales  to  240  m/100  ton^3  as 

indicated  on  Figure  6.5,  and  is  considerably  larger  than  the  mean 

1  /3 

measured  value  of  150  m/100  ton  used  above  in  Equation  (7.1). 

When  the  near-surface  material  does  not  have  a  negligible  tensile 
strength  no  negative  airblast  wing  develops  and  there  is  only  a  coupled 
spall  region.  In  this  instance,  the  radius  of  spall  is  dependent  upon 
the  charge  conf iguration,  and  the  following  equations,  based  on  the  data 
in  Figure  6.5,  can  be  utilized  to  make  predictions: 

HOB:  "R  =  - =  6.5  m/ton1/3  (7.2) 

S  (100)I/3 
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(7.3) 


STS:  Rc  = - ^-77  =  13  m/ton1/3 

S  (lOO)1^ 

STC  or  HBS:  Rc  = - =  19.5  m/to n1/3  (7.4) 

S  (100)i/3 

The  maximum  depth  of  spall  can  be  obtained  from  the  following  equations: 
STS  or  STC:  z$  =  depth  of  water  table*  or  3.2(W,tons)^,m 

whichever  is  smaller  (6.3) 

HBS:  zs  =  depth  of  water  table*  or  5.0(W,tons)^,m 

whichever  is  smaller  (6.4) 

The  high  explosive  predictions  should  be  accurate  to  approximately 
±33  percent. 

C.  Spall  Associated  With  Nuclear  Detonations 

Referring  to  Figure  6.5,  it  can  be  seen  that  R"s  is  approximately 
50m/ (100  ton)^3  for  JANGLE-U  and  140m/ (100  ton)^3  for  the  bermed 
high  explosive  event  presented.  Assuming  that  any  difference  in  the 
spall  data  from  these  two  events  is  attributable  to  the  energy  available 
to  create  ground  shock,  one  obtains  a  nuclear  to  high  explosive  ratio  of 
approximately  36  percent.  R"s  for  the  coupled  nuclear  spall  region  can 
be  estimated  from  Figure  6.5  to  be 

Rs  =  (50) (10)1/3  =  108  m/ktl/3  (7.5) 

In  a  similar  manner  zs  can  be  estimated  from  Equation  (6.3)  as 
z$  =  depth  of  water  table* 

or  (3 .2) (1000) 1/6  =  10.1(W,kt)1/6,m 
whichever  is  smaller.  (7.6) 

*If  there  is  a  rock  layer  at  a  depth  of  zs  or  less,  disregard  the  water 
table. 
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When  the  near-surface  material  has  a  negligible  tensile  strength  a 
negative  airblast  wing  region  can  be  expected.  Definition  of  the  shape 
of  this  region  can  again  be  obtained  through  the  use  of  Figure  7.1.  For 
this  instance,  the  nuclear  "Rs  can  be  estimated  as: 

Rs  =  (0.78) (150) (10)1/3  =  252  m/ktl/3  (7.7) 

Equation  (7.7)  was  obtained  by  taking  the  "R  ratio  of  nuclear  to  high 
explosive  from  Figure  7.1  for  a  constant  depth,  which  is  0.78,  and 
applying  this  ratio  to  Equation  (7.1)  with  appropriate  changes  in  yield. 
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SECTION  VIII 
Summary  and  Conclusions 

This  report  presents  theoretical,  numerical,  and  empirical  analyses  of 
soil  and  rock  spall  due  to  near-surface  explosions,  both  chemical  and 
nuclear. 

Section  VII  presents  an  empirically  based  technique  for  predicting  the 
maximum  radius  of  spall,  R  ,  and  the  maximum  depth  of  spall,  z  ,  for  a 
single  near-surface  explosion,  either  chemical  or  nuclear.  The  prediction 
technique  is  based  on  analyses  of  airblast  and  explosive  ground  motion 
data  presented  in  Sections  V  and  VI. 

Section  III  presents  a  one-dimensional  numerical  analysis  of  spall  in 
a  tensilely  weak,  hysteretic  material  with  a  free  surface,  using  the 
computer  code  STEALTH  ID.  The  results  must  be  considered  preliminary,  but 
they  do  demonstrate  the  ability  of  STEALTH  ID  to  handle  propagation 
problems  involving  the  creation  of  numerous  new  boundaries  (spall  planes), 
as  well  as  the  influence  of  gravity.  One  problem  which  will  arise  in  two- 
dimensional  spall  calculations,  which  does  not  arise  in  one-dimensional 
calculations,  is  how  to  handle  rejoin  grid  mismatch,  i.e. ,  the  fact  that 
two  grid  points  which  coincide  prior  to  spall  probably  will  not  coincide 
after  rejoin.  This  situation  also  arises  in  numerical  analyses  of 
explosive  welding  [Merkle  and  Cannon  (1977)]. 

Section  IV  presents  a  series  of  theoretical  analyses  of  the  pore  air 
effect,  a  principal  near-surface  spall  mechanism  in  airblast  loaded  dry 
soil.  Several  of  these  analyses  yield  finite  difference  equations 
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suitable  for  use  in  parametric  studies.  Even  closed  form  linear  analyses 
demonstrate  the  dramatic  influence  of  both  boundary  conditions  and  soil 
permeability  on  the  transient  pore  air  pressure  distribution  which  can 
cause  spall. 

A  fundamental  understanding  of  soil  spall,  especially  that  caused  by 
local  airblast  loading,  requires  an  effective  stress  approach  in  which 
both  the  soil  skeleton  stress-strain-strength  relations  and  the  transient 
effect  of  flowing  pore  fluid  are  considered  explicitly.  This  approach 
underlies  other  numerical  and  empirical  analyses  which  yield  results 
useful  to  designers  of  hardened  structures  and  their  shock  isolation 
systems. 
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APPENDIX  A 


Isothermal  Fluid  Diffusion  Through  a  Rigid  Porous  Medium 


Darcy's  law  for  percolation  of  a  fluid  through  a  porous  medium  can  be 
written  in  the  form 


3P 

3x 


(A.l) 


where 

v  =  fluid  discharge  velocity  (flow  rate  divided  by  total  area) 

=  effective  permeability 
P  =  fluid  pressure 

x  =  distance 

Under  isothermal  conditions,  the  fluid  is  assumed  to  obey  the  perfect 
gas  law 

P  =  ap  (A. 2) 

where 

a  =  constant 
p  =  fluid  mass  density 

When  the  porosity  of  the  porous  medium  remains  constant,  the  equation 
of  fluid  mass  conservation  takes  the  form 

"  H  *  -  £<»*>  I*'3) 


where 


n  =  porosity,  defined  by  Equation  (E.l) 
Introducing  Equations  (A.l)  and  (A. 2)  into  Equation  (A. 3)  yields 
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•  •-  ■  - 

T  r-.  ■  ^  'V  V"-.'  - - r  . - ; - c - - - - a - - - s - - - - - - . - - - - - ... 

•  :  ■  '  ■  *  -  -  *  • 

(M 

n3P  n  3  ,  PSP  v 
a3t  "  B1  "SPaSx' 

or 

» _!i  a  (P3P} 

3t  ~  n  ^xlK3x; 

(A. 4) 

• 

Equation 

(A. 4) 

can  be  written  in  the  nonlinear  form 

1 

3P  _  B1  322, 

3t  "  2n  .  2l  ; 

3x 

(A. 5) 

• 

On  the  other  hand,  if  the  spatial  variation  of  p  is  ignored  in 

Equation 

(A. 3) 

,  substitution  of  Equations  (A.l)  and  (A. 2)  yields 

where 

3P  B1P  »2P  ^J> 

3t  =  "  9?  =  V 

(A. 6) 

• 

81P 

D"r 

(A. 7) 

Equation 

(A. 6) 

can  be  considered  linear  if  the  parameter  D  is  assumec 

to 

• 

be  constant. 

'  % 

| 

\ 

1 

f » 

[* 

L _ ...  .  . 
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APPENDIX  B 


Transient  Linear  Pore  Air  Diffusion  in  a 

Vented  Layer 

It  is  desired  to  calculate  the  dynamic  pore  air  pressure,  P(x,t) 

f 

satisfying  the  following  conditions: 

3P  n  32P 

‘3r* 0  a? 

(0<x<l) 

(A. 6) 

P(0,t)  =  F ( t) 

( t>0) 

(B.l) 

P(l,t)  =  0 

(t>0) 

(B.2) 

P(x,0)  =  0 

(0<x<l) 

(B.3) 

First,  consider  the  associated  problem  having  homogeneous  boundary 

conditions  [Hildebrand  (1962:431);  Taylor  (1948:229)],  which  is 

3P  32P 

(0< x<2H) 

(B.4) 

P(0,t)  =  0 

(t>0) 

(B.5) 

P(2H, t)  =  0 

(t>0) 

(B.6) 

P(x,0)  =  P^x) 

(0< x<2H) 

(B.7) 

Using  the  method  of  separation  of  variables  [Hildebrand  (1962:430)],  we 
assume  that 

P(x,t)  =  f>(xH(t)  (B.8) 

Substitution  of  Equation  (B.8)  into  Equation  (B.4)  yields 

<njj=  afi"^  (B.9) 

and,  dividing  both  sides  of  Equation  (B.9)  by  the  product  afi\p  =  aP, 
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assumed  not  to  be  zero  everywhere,  we  obtain 

I  i  t>" 

a  <p  ~  T~ 


(B. 10) 


Since  i>  is  a  function  of  x  only,  and  »);  is  a  function  of  t  only,  and  yet 
Equation  (B.10)  must  hold  for  all  values  of  x  and  t,  it  must  be  that  both 


sides  of  Equation  (B.10)  are  equal  to  the  same  constant,  i.e., 

1  jL  t"  2 
a  IT 


(B. 11) 


Equation  (B.ll)  yields 


+  Jt  =  0 


(B. 12) 


+  aa>  i|>  =  0 


(B. 13) 


The  solutions  to  Equations  (B.12)  and  (B.13)  are 

£(x)  =  cos  ux  +  C2  sin  u> 


(B.14) 


4>(  t )  - 


-au2t 


(B  .15) 


so  that  Equation  (B.8)  yields 


P ( x , t)  =  0(x)^(t)  =  (C^cos  u>x  +  C^sin  u>x)e 


2 

-aoi  t 


(B  .16) 


where 


C4  =  C1C3 


(B.17) 


C5  =  C2C3 

In  order  to  satisfy  Equation  (B.5)  it  must  be  that 

c4  =  0 


(B. 18) 


(B.19) 
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and  in  order  to  satisfy  Equation  (B.6)  it  must  be  that 

2wH  =  nir  (n  =  1,2,3,  ...) 


( B .  20 ) 


which  means  that 


nir 

"  =  ?H 


(n  —  1,2,3,  » . . ) 


(B.21) 


Thus  an  infinite  number  of  functions  of  the  form 


2  2 

n  *  at 


r  nirX  . 

C5  Sin  “7H  e 


each  satisfies  Equations  (B.4),  (B.5)  and  (B.6),  and  since  Equation  (B.4) 
is  linear  and  Equations  (B.5)  and  (B.6)  are  homogeneous,  any  linear 
combination  of  the  above  functions  is  also  a  solution.  Thus  we  can  write 


uu 

:,t)  =  E 


2  2 

n  H  at 


Bn  sinJ3fHe 


(B.22) 


The  constants  Bn  (n  =  1,2,3,  ...)  are  determined  from  Equation  (B.7). 


00 

P(x,0)  =  P.(x)  =  E  Bn  sin 


(B  .23) 


Equation  (B.23)  is  a  Fourier  sine  series,  and  since,  with 


(B.24) 


we  have 


,  m*x  nirX  . 

f  sin  — ^  sin  — dx 


=  —  f  sin  me  sin  ne  de 


=  -  /  [cos(m-n)e  -  cos(m+n)e]de 
n  0 
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and 


2/  sin2  S  dx  =  —  /  sin2ne  do 
0  w  0 


—  /  (1  -  cos  2ne)de  =  H 
11  0 


(B.26) 


then  if  both  sides  of  Equation  (B.23)  are  multiplied  by  sin  rmrx/2H,  and 
the  results  integrated  on  x  from  0  to  2H,  the  result  is 


s  Pl(x)  sin^dx  =  BmH 


so  that 


Bn  -  ff  l  Pi<*>  s1n  Tff 


(B.27) 


We  now  consider  the  problem  identical  to  that  defined  by 
Equations  (A. 6),  (B.l),  (B.2)  and  (B.3),  except  that  the  function  F(t)  in 
Equation  (B.l)  is  assumed  to  be  constant,  i.e., 

P(0,t)  =  F  (t>0)  (B .28) 

The  approach  to  this  problem  is  to  assume  that  [Carslaw  and  Jaeger 
(1959:99)] 


P  =  U  +  W 


(B.29) 


where 


4.0 

dx^ 

(0<x<l) 

(B.30) 

U(0,t)  =  F 

(t>0) 

(B.31) 

U(l,t)  =  0 

(t>0) 

(B.32) 
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and 


3W  n  a  w 

3t  '  3X2 

(0<x<  1 ) 

(B.33) 

w(0,t)  =  0 

(t>0) 

(B.34) 

w(l.t)  =  0 

(t>0) 

(B.35) 

W(x,0)  =  — U ( x) 

( 0<  x<  1 ) 

(B.36) 

is  easily  found 

to  be 

U(x)  =  F (1  - 

T> 

(B.37) 

and  the  expression  for  W  is  found  from  Equation  (B.22) 


.»  )  n  nitx 

:,t)  =  L  Bn  sin  — r- 


n2ir2Dt 


(B.38) 


where  Equations  (8.27),  (B.36)  and  (B.37)  yield 


D  2F  .  n  Xx  .  nnx  , 

Bn  *  -~T^  (1  -T)  ""  T  dl 


(B.39) 


The  integral  on  the  RHS  of  Equation  (B.39)  can  be  integrated  by  parts  by 
setti ng 


u  =  1  - 


a  dx 

du  y 


dv  =  sin  dx 


,,  1  , _ nirx 

V  =  -  —  COS  — T- 
nn  1 


(B.40) 

(B.41) 

(B.42) 

(B.43) 


so  that 


APPLIED  RE/EAPCH  A//OCIATE/finC. 


(B.44) 


IB  ,  1-1 

_ n  In  x,  nwx  |  1  ,  nirx  . 

-  -7F  -  -  m(1  -  7>  cos  T  I  -n—  £  cos  T  dx 

_L 

~  nir 

and  therefore 

B  .-it 
n  nir 


(B.45) 


Thus,  Equations  (B.29),  (B.37),  (B.38)  and  (B.45)  yield 


P(x,t) 


=  F 


oo 


2  2 
n  »  Dt 


(B.46) 


[cf.  Carslaw  and  Jaeger  (1959:103)]. 

We  can  now  return  to  the  original  problem  defined  by  Equations  (A. 6), 
(B.l),  (B .2)  and  (B.3),  for  which  the  boundary  value  in  Equation  (B.l)  is 
a  function  of  time.  Because  the  system  is  linear,  the  solution  can  be 
written  as  a  Duhamel  integral  in  the  form 


t 

P(x,t)  =  /  F(x)h(x,t-x)dx  (B.47) 

0 


where 

h(x,t-x)  =  solution  for  a  unit  impulse  at  (0,x). 

However,  what  we  have  obtained  above  is  not  h(x,t-x),  but  its  integral, 

•  s(x,t-x)  =  solution  for  a  unit  step  at  (0,x). 

In  order  to  use  the  step  function  response,  we  integrate  Equation  (B.47) 
by  parts,  by  setting 

t  u  =  F(x)  (B.48) 

du  =  F'(x)dx  (B .49) 
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dv  =  h(x,t-x)dx 
v  =  -s(x,t-x) 


(B.50) 
(B  .51) 


so  that 


t  t 

P(x,t)  =  -F(x)s(x,t-x)  |  +  f  F' (x)s(x,t-x)dx 

0  0 


=  -F(t)s(x,0)  +  F(0)s(x,t)  +  /  F' (x)s(x,t-x)dx 

0 


(B .52) 


where,  from  Equation  (B.46),  we  have 


s(x,t-x)  = 


(1  -4) 


2  y  i  . 

-:L  -  si 


n27r2D(t-xJ 


*  n=l  n 


_  -  nitx  _ 
sin  — i—  e 


(B.53) 


s(x,t)  -  (1  -4) 


uu 

2  y  1  .  nu; 

- —  sin  — y 

ir  _  i  n  1 


2  2 

nVDt 


(B. 54) 


Now  if 


so  that 


F(x)  =  P0e“aA 


(B.55) 


F  ( 0 )  =  P, 


(B. 56) 


F'U)  -  -aP0e" 
then  the  terms  in  Equation  (B.52)  are 


F (0)s(x,t)  =  PQ  (1 


GO 

x)  2  Z  1 

■  v  ~  *  r:  n 


_  n2ir2Pt 
sini^e  ,2 


(B.57) 


(B.58) 
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t  t 

/  F'(x)s(x,t-x)dx  =  -aP  (1  -  4)  /  e-aX  dx 
0  0  1  0 


-  p0d  -  ■y)(e"at  - 1) 


(B .59) 


Finally,  substitution  of  Equations  (B.58)  and  (B.59)  into  Equation  (B.52) 
yields  r 
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If  we  set 


r® 


Dt 


=  T 


l2a 


=  8 


4  =  5 


then  Equation  (B.60)  can  be  written  in  the  form 


(B.61) 

(B.62) 
( B . 63 ) 


P(C.T)  =  P, 


(1-Oe 


-ST 


- 1 T.  i  =,• 

’  n=l  " 


.  _  ,  n2,2e-"2’2T  -  Be~eT 

sin  nir£  (  - - 

n  Sr  -  8 


(B.64) 
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APPENDIX  C 


Transient  Linear  Pore  Air  Diffusion  in  an  Unvented  Layer 

It  is  desired  to  calculate  the  dynamic  pore  air  pressure,  P(x,t), 
satisfying  the  following  conditions: 


3P  n  32P 

^  =  V 

(0<x<l) 

(A. 6) 

P(0,t)  =  F(t) 

(t>0) 

(c.l) 

Hd^)  -° 

(t>0) 

(C.2) 

P(x,0)  =  0 

( 0<  x<  1 ) 

(C.3) 

The  above  solution  can  be  most  easily  obtained  as  the  first  half  of  the 
solution,  symmetric  about  x  =  1,  for  which 


3P  n  3  P 

(0<x<21) 

(C.4) 

P(0,t)  =  F(t) 

(t>0) 

(C.l) 

P(21,t)  =  F  ( t) 

(t>0) 

(C.5) 

P(x,0)  =  0 

(0<x<21 ) 

(C.6) 

We  will  again  use  the  solution  for  homogeneous  boundary  conditions 
defined  by  Equations  (B.22)  and  ( B.27 ). 

We  now  consider  the  problem  identical  to  that  defined  by 
Equations  (C.4),  (C.l),  (C . 5 )  and  (C . 6) ,  except  that  the  function  F(t)  in 
Equations  (C.l)  and  (C.5)  is  assumed  to  be  constant,  i.e. , 

P(0,t)  =  P(21 , t)  =  F  (t>0)  (C.7) 
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As  in  Appendix  B,  we  assume  a  solution  of  the  form 


P  =  u  +  w 

(C.8) 

where,  in  this  case 

4=0 

dx^ 

(0<x<21 ) 

(C.9) 

U(0.t)  =  F 

(t>0) 

(C.10) 

U(21, t)  =  F 

(t>0) 

(C.ll) 

and 

9W  n  92W 
^  3x2 

( 0<  x<  2 1 ) 

(C.12) 

W(0,t)  =  0 

(t>0) 

(C.13) 

W(21, t)  =  0 

(t>0) 

(C.14) 

W(x,0)  =  -U(x) 

(0<x<21 ) 

(C.15) 

The  expression  for  U  is  easily  found  to  be 

U(x)  =  F  (C.16) 

and  the  expression  for  W  is  obtained  from  Equation  (B.22) 

°°  _  n2»2Pt 

W(x,t)  =  Y,  Bnsin-^Jje  41?  (C.17) 

n=l 

where  Equations  (B.27),  (C.15)  and  (C.16)  yield 
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cf.  [Taylor  (1948:232)].  Thus  Equations  (C.8),  (C. 16) ,  (C.17)  and  (C.18) 
yield 


P(x,t)  =  F 


1  - 


il 

11  n=l 


(2n-l)irX 
sin - — 

— 


_  (2n-l)2ir2Dt 
.  412 


(C.19) 


We  can  now  return  to  the  original  problem  defined  by  Equations  (C.4), 
(C.l),  (C.5)  and  (C.6),  for  which  the  boundary  value  in  Equations  (C.l) 
and  (C.5)  is  a  function  of  time.  Equation  (B.52)  again  applies,  where 
from  Equation  (C.19), 


s(x,t-x)  =  1  - 


il 

11  n=l 


.  (2n-l)irx 

sin - _1 - 

— am. 


(2n-l)2Tt2D(t-x) 

412 


(C.2 0) 


2  2r 


s(x,t)  =  1  - 


iZ 

*  n=l 


sin 


(2n-l)trx 


21 

2n  -  T 


41 


(C.21) 


Assuming  the  boundary  value  inputs  are  again  defined  by  Equations  (B.55), 
(B .56)  and  (B.57),  the  ter.r.s  in  Equation  (B.52)  are 


and 


2r 


F(0)s(x,t)  =  Pr 


f  s1„(2n-D.x  -  (2n-1)  Dt  1 

,  4  >  Sin - 2 T~  „  412 

1  '  '  n.l  11  8 


(C.22) 


f  F’(x)s(x,t-x)dx  =  -oPn  /  e  aX  dx 
0  0  0 


„  ”  .  (Zn-lUx  - 

4aP  y  sin  ■'  41^  t 

*  —  ^  -  e  1 


(2n-l)2ff2D 


-  a 


41 


n=l 
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r  t 

=  p«  e~aX 

0  L  Jo 


4ap„  V 
+  — *  L 


.  (2n-l)*x 

>  in - ?5t — 


(2n-l)Vpt 


(2n-l)V 

:  41* 


-  o  IX 


11  n=l  ~ 


(2n-l)VD 


=  p0(e  at  -  1} 


ia 


( 2n— 1 )  ttx 


(2n-l)27r2Dt 


(2n-l)2iT2D 


(C.23) 


Finally,  substitution  of  Equations  (C.22)  and  (C.23)  into  Equation  (B.52) 
yields 


p(x,t)  =  Poe- 


4Pq  y  sin 
~V  2n  -  1 


(2n-l)2ir2D 


'2n-l)2ir2Dt 


41  -at 

-ae 


(2n-l)  t  D 


(C  .24) 


cf.  [Carslaw  and  Jaeger  (1959:105)], 
If  we  set 


(C.25) 


(C .26) 
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then  Equation  (C. 24)  can  be  written  in  the  form 


PU.T)  =  PQ  e  BT 


sin  Nir£ 


V,2e-N2’2T  -  8e-6T  Y| 

,  »2,2  -  ,  ) 


(C.29) 
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APPENDIX  D 


Adiabatic  Fluid  Diffusion  Through  a  Rigid  Porous  Medium 


Darcy's  law  for  percolation  of  a  fluid  through  a  porous  medium  can 
written  in  the  form 


v  =  -B,— 
lax 


(A.l) 


where  the  terms  in  the  above  equation  are  defined  in  Appendix  A. 

Under  adiabatic  conditions,  the  fluid  is  assumed  to  obey  the  relation 

P  =  aPY  (D.l) 

where  y  =  1.4  for  air,  and  the  remaining  terms  are  defined  in 
Appendix  A. 

When  the  porosity  of  the  porous  medium  remains  constant,  the  equation 
of  fluid  mass  conservation  takes  the  form 


it 


(A. 3) 


where  n  is  porosity,  as  defined  by  Equation  (E.l).  Introducing  Equations 
(A.l)  and  (D.l)  into  Equation  (A. 3)  yields 

n/j\Y  I  _  r  JL  [r£.^Y  9P 
y'  &  *  a  at  -  Biax  V  ax 


ap  Vp 
at  = 


1-1/y 


-  ( 
ax  l 


,i/y  sp 


Setting  y  =  1  in  Equation  (D.2)  yields 


(D.2) 
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<• 


dP  j_  /pjPs 

9t  n  gx  '  sx' 


(A.4) 


and  neglecting  the  spatial  variation  of  p  (i.e.,  P^Y)  in  Equation 
( D.2)  yields 


9P  _ 

9t 


biYp 


1  T‘  32P 

7T~? 

3x 


(D.3) 


where 

Bi  YP 

D  =  —p  (D.4) 


Equation  D.3)  can  be  considered  linear  if  the  parameter  D  is  assumed  to 
be  constant. 

If  no  assumption  is  made  concerning  the  terms  in  Equation  (D.2),  and 
the  RHS  is  expanded,  the  result  is 


9P 

3t 


biYp 


1-1/y  r 


I  py 
r 


I-  1 


,dP)2  +  pl/Y 
^  9X2 


1 

n 


?  2  ' 
/9P>  +  p9  P 

(^’  ,P^7 


(D.5) 


which  is  the  nonlinear  equation  used  by  [Zernow,  et  al  (1973:101)].  The 
first  term  in  brackets  on  the  RHS  of  Equation  (D.5)  reflects  the  spatial 
variation  of  p;  the  second  term  reflects  the  spatial  variation  of  v. 

S;_e  Appendix  F  for  a  finite  difference  approximation  to  Equation 
(D.5). 
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APPENDIX  E 

Propagation  Velocity  in  Spalled  Soil 

It  is  assumed  that  the  soil  skeleton  is  distended,  and  therefore 
occupies  volume  but  carries  no  load.  Thus  the  material  can  transmit 
hydrostatic  pressure,  but  not  shear. 

The  standard  soil  phase  diagram  applies. 


and  the  following  definitions  apply: 

porosity: 

total  density: 

volumetric  strain: 


(E.l) 


(E.2) 


(E.3) 
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The  sound  speed,  c,  at  pressure,  P,  for  the  mixture  of  soil  particles 
in  air  is  give  j  the  equation 


r2  1  dP 
c 


(E.4) 


where 


VT  -  VT0(1  +  «) 


(E.5) 


dVT  =  VT0de 


(E.6) 


PT  -  VT0(l+e)  = 


(E.7) 


Thus,  Equation  (E.4)  can  be  written  in  the  form 


_2  (1  +  E  r  dp 

c  = 


(E.8) 


For  rapid  pressure  changes,  the  relation  between  pore  air  pressure  and 
pore  air  (void)  volume  is  the  adiabatic  equation 


p-p°ft)W¥ 


(E.9) 


where  y  =  1.4  for  air.  Therefore,  we  have 


o  o 


n  +e\  (Y+1^ 


n„  n„ 
o  \  o 


(E.10) 


so  that  Equation  (E.8)  can  be  written  in  the  form 


2  yPo  flJ.  v2  (  no 


(E.ll) 
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and  finally 


(E.12) 


The  value  of  c  for  a  typical  soil  has  been  tabulated  as  a  function  of 
strain  by  [Merkle  (1980:38)]. 
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APPENDIX  F 


Three-Dimensional  Adiabatic  Fluid  Diffusion  Through  a  Rigid,  Porous, 
Isotropic  Medium  "  " 

Darcy's  law  for  fluid  diffusion  through  a  porous,  isotropic  medium  can 
be  written  in  the  form 


v  = 


-B1VP 


(F.l) 


where 

v 


P 

V 


=  fluid  discharge  velocity  vector 
=  effective  isotropic  permeability 
=  fluid  pressure 


=  gradient  vector  operator 


xl’  x2’  x3  =  rectangular  Cartesian  coordinates 


ep  e2»  e^  =  unit  vectors  in  the  three  coordinate  directions 


Under  adiabatic  conditions,  the  fluid  is  assumed  to  obey  the  relation 

P  =  aPY  (F.2) 

where 

a  =  constant 

p  =  fluid  mass  density 

y  =  ratio  of  specific  heats  (=  1.4  for  air) 

When  the  porosity  of  the  porous  medium  remains  constant,  the  equation 
of  fluid  mass  conservation  takes  the  form 
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01 


n||.  -V  •  (pv) 


(F.3) 


where 


n  =  porosity 


Introducing  Equations  (F.l)  and  (F.2)  into  Equation  (F.3)  yields 


i-i 


3P  Vpl~1/T  V  ■  (pl'r  VP) 
3t  n 


(F.4) 


Setting  y  =  1  in  Equation  (F.4)  yields 


||  .  4  v-(PVP) 


(F.5) 


and  neglecting  the  spatial  variation  of  p  (i.e.,  in  Equation  F.4 


yields 


aP  ?  2 

3t  n 


(F.6) 


where 


(F.7) 


Equation  [F.6]  can  be  considered  linear  if  the  parameter  0  is  assumed  to 
be  constant. 

If  no  assumption  is  made  concerning  the  terms  in  Equation  (F.4),  and 
the  RHS  is  expanded,  the  result  is  [Hildebrand  (1962:278)] 
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BlYP 


1-1 /y 


+  pl/^y^p 


(F.8) 


The  first  term  in  parentheses  on  the  RHS  of  Equation  (F.8)  reflects  the 
spatial  variation  of  p;  the  second  term  reflects  the  spatial  variation  of 
"v. 


It'  we  set 


piJ,k  -  P^ure  at  (x,,^) 

then  the  two  dimensional  finite  difference  approximation  to  Equation  (F.8) 
is  [Crandall  (1956:246,376)] 


Pi,j,k+1  Pi,j,k  _  B1 


At 


(Pi+l,j,k  ~  Pi-l,j>k)2+  (Pi,j+l,k  ~  Pi,j-l,k) 


4h 


+  Yi,j,k 


Pi-l,j,k  *  Pi+l,j,k  +  Pi,j-l>k  +  Pi,j+l,k  4Pi,j,k 


(F.9) 


where 


h  =  ax  =  Ay 


(F.10) 


so  that 


B,  At 

p  —  p  +  ^ 

i,j,k+l  i,j,k 


nh 


^Pi+l,j,k  “  Pi-l,j,k^  +  ^Pi,j+l,k  Pi,j-l,k^ 
-  ? - 


+  YPi,j,k  ^Pi-l,j,k  +  Pi+l,j,k  +  Pi,j-l,k  +  Pi,j+l,k  4Pi,j,k 


(F .11) 


For  the  one  dimensional  case.  Equation  (F.ll)  reduces  to 
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<• 


B,  At 

k+1  =  k  +  7 

1,K  1  ’  nlr 


(-i,'V-i,k)^<,lt(l’m,k  -  2Pi,k *  pi-i,k)‘ 


(F.12) 


An  alternate  formulation  of  the  above  equations  results  in  a  pseudo- 
linear  partial  differential  equation  without  any  approximation. 
Substituting  Equation  (F.l)  into  Equation  (F.3)  yields 


n  If  =  -v  •  (— pBjVP )  =  B1V  •  (pVP) 


Now  Equation  (F.2)  can  be  written  in  the  form 


p  - 


so  that  substituting  Equation  (F.14)  into  Equation  (F.13)  yields 


_3_ 

3t 


=  BjV  * 


(-)1/tvp 


or 


(p1/y)  V  •  (P1/tVP) 


Now  on  the  RHS  of  Equation  (F.15)  we  have 


p1/tvp  =  r/---V-i  v(p1/y+1) 

so  that  Equation  (F.15)  can  be  written  in  the  pseudo-linear  form 


The  one-dimensional  form  of  Equation  (F.17)  is 


(F.13) 


(F.14) 


(F.15) 


(F.16) 


(F.17) 
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£(e1/v)  - 


3X 


(F - 18) 


For  thr  isothermal  case,  setting  Y  =  1  reduces  Equation  (F.18)  to 
Equation  (A. 5).  The  finite  difference  approximation  to  Equation  (F.18)  is 


p 1/y  p1/y 

*i,k+l  _  *i,k 


n  1/y  + 


1+1  l+i  l+i 

pY  _  2dY  +  pY 

1-1, k  f  i,k  Ki+l,k  (F.19) 


so  that 


1/v  1/v  MAt)  ^  +  1  ±  +  1  i  +  1 

ik+1  =  Pik  +  - - - - 2  (Pi-l  k  "  2Pi  k  +  Pi+i  k)  (F*2°) 

1  1,K  n  1/y  +  IHaxT  1  i,k  1,k  1+1»k 


1+1  1  + 


'(1/y  +  1)  (ax) 


where 


pi/ Y+1  _  /pl/Y  4+Y 
Ki ,k+l  ~  'Ki,k+l' 


(F.21) 


p  _  (p!/Y  \Y 

v i , k+1  "  ^i,k+lj 


(F.22) 


A  similar  pseudo-linear  formulation  can  be  constructed  with  p  as  the 
dependent  variable.  Substituting  Equation  (F.2)  into  Equation  (F.18) 


yields 


n  vWr*1] 


_9p  _ 

at  ~  7R 


(F.23) 


I 
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The  one-dimensional  form  of  Equation  (F.23)  is 


3x 


(F.24) 


The  finite  difference  approximation  to  Equation  (F.24)  is  straightforward, 
but  not  quite  as  convenient  for  computation  as  Equations  (F.20),  (F.21) 
and  (F.22),  since  pressure  boundary  conditions  must  be  converted  to 
density.  Equations  (F.20),  (F.21)  and  (F.22)  are  even  convenient  for  hand 
computation. 
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APPENDIX  G 


Transient  Linear  Pore  Air  Diffusion  in  a  Vented  Layer  Subjected  to  Surface 


Airblast  Loading  Having  a  Negative  Phase 


All  the  equations  of  Appendix  B  apply,  up  to  and  including  Equation 


(B.54).  In  place  of  Equation  (B.55)  we  assume  that 


F(x)  =  P  e-otX(l  -  -pO 

Ao 


(G.l) 


so  that 


F(O)  =  P, 


(G.2) 


F ‘ ( x )  =  -«P0e  aX(l  -  -£-)  e' 

o  o 


=  -a ( 1  +  -4-)Pne~aX  '  IT  PnXe_aX 
aXo  0  X0  ° 


(G .  3 ) 


If  we  set 


FJ(X)  .  -.(1  -  ^)P0e-“A 


(G.4) 


F.(x)  =^P0xe"“x 


(G .  5 ) 


then  Equation  (G.3)  can  be  written  in  the  form 


F ’ ( x )  =  Fj(x)  +  F’(x) 


(G.6) 


The  terms  in  Equation  (B.52)  then  are 
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F(0)s(x,t)  = 


e0  (l-^) 


2  2 

n  it  Dt 


*  n=l  n 


sin  e 


(6.7) 


/  Fj( x) s(x, t-x)dx  =  PQ(1  +  ^-)(1  -  ^ 
0  o 


+  ^-)d  -  4)(e_at  - 


2(“  +  ^)P°  f.  1  • 


2  2 

-  n St  Dt 


T  -at 

-  e 


2  2 
n  *  D 


(G.8) 


and  finally 


/  Fl(x)s(x,t-x)dx  = 

c\ 


po  TT(1  -  T>  S  ^  di 

0  Ao  1  0 


E  i  si, 


2  2 

n  V  Dt 


0  )  i  nirx  „ 
.  /  i  sin  i  e 

tt  x  i  n  1 
o  n=l 


-2  t  \ 
/  xe 
0 


2  2 

nVD 


(G.9) 


Both  integrals  in  Equation  (G.9)  are  of  the  form 


I  =  /  xerxdx 


(G. 10) 


which  can  be  integrated  by  parts  by  setting 


U  =  x 


(G.ll) 


dU  =  dx 
dV  =  erx  dx 


(G.12) 
(G . 13) 


V  =  -p  erX 


(G.14) 
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AD-A125  592 


UNCLASSIFIED 


EFFECT  OF  SPALL  ON  THE  CHARACTERISTICS  OF  EXPLOSIVE  ' 
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so  that 


I  =  ^  erx|  -  j  f  erx  dx 
0  r  0 


/ a  1»  rx 

<7  -  7)e 


1  -  (7  -  4)ert  +  4 

A  1  w.L  pU 


(G.15) 


Thus,  Equation  (G.9)  yields 


r  F^(»)s(x,t-»)dx  =  p0  4  -  <1  *  4)e 


1  \  -at 


00 

(  n2ii2Dt 

"  l2 
e 

2“Po 

y  1 

•  n  nirx  . 
sin  --j—  < 

t 

1 

e~at 

wXo 

h" 

/  n2ir20  \2 

IV  i2  ) 

n2tr2D 

,2  -  1 

( n2ii2D  \2 

V  l2  /  J 

j 

When  Equations  (G.7),  (G.8)  and  (G.16)  are  substituted  into  Equation 


(B.52),  the  coefficients  of  key  functions  are 


P0(l  -7):  1  -  e 


-  1  + 


.-at 


aA, 


aX„  aA„ 
0  O 


and  with 


„Vd 


r 


n 


2P 

IT 


(G - 17) 


(G. 18) 
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rn  -  “  +  *0(rn  - 


(G.19) 


2P  °  T~ 

0  1  nirX  „-at .  0  a  j.  at 

~r<  ' 


_n _  +  at 


rn  “  x0(rn  -  a)2  xo{rn 


(G.20) 


Thus,  if  we  set 


_ 9L_  +  _ D _ _  _  e 

r„  -  a  ~TZ  xT  n 
n  x0(r„  -  a) 


(6.21) 


then  Equation  (B.52)  takes  the  form 

P(x,t)  =  PQ  f  (1  -4)(1  -  ^)< 

L  o 


-  -  X!  4  sin 


n  n=l  n 


'  -r  t 

^  (i  +  «n)e  " 


- 

at  "|  -at 

e  J 


(G.22) 


If  we  set 


(B.61) 


(B.62) 


(G.23) 


T-« 


(B.63) 
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then 


APPENDIX  H 


The  volume  change  of  the  soil  skeleton  depends  on  effective  stress. 

p  -  u  «  k$x  (H.l) 

The  volume  change  of  the  initial  volume  of  pore  fluid  depends  on  the  pore 
pressure. 

u  =  k^  (H.2) 

The  rate  at  which  new  pore  fluid  flows  into  the  soil  void  space  depends  on 
the  pore  pressure  difference. 

(y  -  X)  =  (p  -  u)  (H.3) 

Written  together.  Equations  (H.l),  H.2)  and  (H.3)  take  the  form 
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(M 


ft 


ft 


V 


-Cvx  +  Cvy 


kgX  +  u  =  p 

-  u  =  0 
+  u  =  p 

Eliminating  the  pore  pressure,  u,  we  obtain 

Cvx  -  Cvy  +  ksx  =  0 

-Cvx  +  Cvy  +  kwy  =  p 

or,  in  matrix  form 


(H.4) 


(H.5) 


ks 

0 


kw 


1  -1 
L-1  l 

In  operator  notation.  Equations  (H.5)  take  the  form 

(CVD  +  ks)x  -  CvDy  =  0 
-Cv0x  +  (CVD  +  kw)y  =  p 

We  eliminate  y  as  follows: 


P 

v  / 


(H.6) 


(H.7) 


(CVD  +  ks)(CyD  +  kw)x  -  CvD(CyD  +  kw)y  =  0 


(+) 


-(CyD)4x 


+  CvD(CvD  +  kw)y  =  CvDp 


[Cv(ks  "  kw)D  +  k$kw]x  =  CyDp 


(H .  8) 


Wher  Equation  (H.8)  has  been  solved  for  x,  the  first  of  Equations  (H.4) 
yields  the  solution  for  u: 


u  =  p  -  k$x 


and  the  second  of  Equations  (H.4)  yields  the  solution  for  y; 

i 

k. 


y  =  7T 


(H.9) 


(H.10) 


w 


Equation  (H.8)  can  be  written  in  the  form 


D  + 


kskw 

♦  y 


X  =  ( 


k  +  k 
KS  Kw 


)Dp 


(H.ll) 


I 
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or 

(D  +  a)x  =  bDp  (H.12) 

where 


■ 

a  ■  +_V 

(H.13) 

• 

6  ‘  k<  +  kw 

s  w 

(H.14) 

• 

The  complete  solution  to  Equation  (H.12)  for  which  x(0)  =  0 

is  [Cheng 

(1959:15)] 

x(t)  =  b  /  ^  e_a(t_x)  dx  (t>0) 

(H.15) 

• 

0  A 

The  integral  on  the  RHS  of  Equation  (H.15)  can  be  integrated 

by  parts  by 

setting 

• 

U  .  e'a(t-j) 

(H.16) 

dU  =  ae'a*t_x)  dx 

(H.17) 

• 

dV.^dx 

(H.18) 

V  =  p(x) 

(H.19) 

<• 

so  that 

x(t)  =  bp(x)e_a(t~x)  \  -  ab  /  p(x)e"a(t~x)  dx 

0  0 

=  bp(t)  -  bp(0)e_at  -  ab  /  p(x)e“a^t-x)  dx 

0 

(H.20) 

• 

When 

O 

n 

O 

CL 

(H.21) 
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■"  r 


Equation  (H.20)  reduces  to 


x(t)  =  b  p(t)  -  a  /  p(x)e 

0 

When  p(t)  is  a  step  pulse,  i.e.,  when 


dx 


(t>0) 


u(t)  =  P 


1  - 


-at 


T  +~F 
*s  \ 


(t>0) 


and  Equation  (H.10)  yields 


y(t)  =  TT- 


1  - 


s  \  -at 

k  nr  ' e 

s  w 


(t>0) 


Note  that  Equations  (H.24)  and  (H.26)  yield 


x(0+)  =  y(CH-)  =  y~T  T~ 
s  Kw 


(H.22) 


• 

p(t)  *  0 

(t<0) 

=  p 

(t>0) 

(H.23) 

lit =  p6(t) 

• 

then  Equation  (H.15)  yields 

• 

x(t)  =  bP  /  «(x)e"a(t“x^  dx 

0 

un„-at  P  „-at 

“  bPe  -  F~+~k"  6 

Ks  Kw 

(t>0) 

(H.24) 

• 

Equation  (H.9)  then  yields 

(H.25) 


(H.26) 


(H.27) 
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(# 


# 


The  reasonableness  of  the  above  solution  can  be  examined  by 
constructing  the  model  associated  with  Equations  (H.5)  or  (H.6),  which  is 
shown  below. 


Initial  application  of  the  step  pulse  P  produces  equal  instantaneous 

displacements  x(0+)  and  y(0+),  so  that  no  dashpot  force  acts.  The 

springs  k  and  k  resist  the  load  in  parallel.  At  t  =  °°  the  spring 
I  s  w 

k  resists  the  entire  load,  and  the  soil  piston  has  returned  to  its 
w 

original  position. 

^  The  response  is  composite  compression  followed  by  recovery.  Initial 

undrained  loading  produces  a  sharing  of  internal  stress  between  the  soil 
skeleton  and  the  pore  fluid,  in  such  a  manner  that  the  initial  volume 
^  changes  of  soil  skeleton  and  pore  fluid  are  equal.  Subsequently,  drainage 

occurs  in  response  to  the  difference  between  the  back  pressure  and  the 
pore  pressure,  as  a  result  of  which  the  pore  pressure  eventually  becomes 
equal  to  the  back  pressure.  At  this  point  the  change  in  effective  stress 
in  the  soil  skeleton  is  zero,  and  its  initial  undrained  volume  change  has 
been  recovered. 

* 


t 
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When  p(t)  is  an  impulse,  i.e.. 


t 


p(t)  =  6(t)  (H.28) 

then  Equation  (H.22)  yields 

h(t)  =  -abe‘at  (H.29) 

The  fact  that  an  upward  displacement  results  from  a  downward  impulse  can 
be  explained  by  the  fact  that  a  downward  impulse  can  be  viewed  as  the 
successive  application  of  two  step  pulses,  the  first  downward  and  the 
second  upward. 

At 


f(t) 


Equation  (H.24)  therefore  yields 

h(t)  =  bP[e-at  -  e"a(t_At)] 

-  bP[^|  (e-at)at] 

=  -abe"at(PAt)  =  -abe_at  (H.30) 
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The  net  upward  displacement  following  a  downward  impulse  is  due  to  the 
fact  that  displacement  recovery  begins  immediately  following  the 
application  of  a  step  pulse.  This  is  caused  by  flow  of  pore  fluid  into 
the  soil  skeleton  in  response  to  the  difference  between  the  backpressure 
and  the  pore  pressure. 
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APPENDIX  I 


A  =  1 


The  equation  of  motion  for  the  soil  skeleton  is 

{p  -  u)  -  Ksx  =  Msx  (1.1) 

The  equation  of  motion  for  the  pore  fluid  is 

u-Kwy=Mwy  (1.2) 

The  rate  at  which  new  pore  fluid  flows  into  the  soil  void  space  depends  on 
the  pore  pressure  difference 

3Y  (y  -  x)  =  ^  (p  -  u)  (1.3) 

Written  together.  Equations  (1.1),  (1.2)  and  (1.3)  take  the  form 
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+  Kgx 


+  u  =  p 


+  IV 


-  u  =  0 


(1.4) 


-V  +  Cvy 


u  =  p 


Eliminating  the  pore  pressure,  u,  we  obtain 


M$x  +  Cyx  -  C vy  +  K$x 


V  -  Cv*  +  Cvy  +  V  =  P 


(1.5) 


or,  in  matrix  form 


0  x 


0  Mw  y 


1  -l  x  Ks  Ox 

+  cv 

-1  1  y  o  Kw  y 


In  operator  notation,  Equations  (1.5)  take  the  form 

Lix  -  L2Y  =  0 


where 


-L2X  +  L3y  =  p 


Li  =  MSD2  +  CVD  +  Ks 


L2  =  CVD 

L3  =  MWD2  +  CVD  +  Kw 

To  eliminate  y  from  Equations  (1.7),  we  proceed  as  follows: 

L1L3x  “  L2L3y  =  0 
(  +  )  "  L2X  +  l2l3y  =  L2P 

(1^3  -  L^)x  =  L2p 


(1.6) 


(1.7) 


(1.8) 

(1.9) 

(1. 10) 


(1. 11) 


Expanding  Equation  (1. 11),  we  obtain 
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[MsD2  +  CVD  +  K$){MwD2  +  CyD  +  Kw)  -  C2D2]x  =  CyDp 


or 


([MsMw°4  +  CV(MS  +  Mw)D3  ^  (M$Kw  +  M^JD2 


+  VKs  +  Kw>D  +  KSKw>  "  CvDp 


(1.12) 


When  Equation  (1.12)  has  been  solved  for  x,  the  first  of  Equations  (1.4) 
yields  the  solution  for  u: 


u  .  p  -  (M$D£  +  Ks)x 


and  the  third  of  Equations  (1.4)  yields  the  solution  for  y: 


i  t 

y(t)  =  x(t)  +  A  /  [p(x)  -  u(x)]dx 

rs 


'v  0 


Equation  (1.12)  can  be  written  in  the  form 


D4  + 


or 


MM 
S  W' 


(D4  +  bD3  +  cD2  +  dD  +  f)x  =  gDp 


I  M$Mw  * 


/KsKw 

D  +  (  )>x  = 


,MsV 


Dp 


where 


b  = 


CV(M$  -  Mw) 


(1.13) 


(1.14) 


(1.15) 


(1.16) 


(1.17) 
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# 


c 


MsKw  +  MwKs 
MsMw 


d  = 


CV(KS  +  V 

MsMw 


f  = 


K,K 
s  w 

MM 
s  w 


9  = 


MM 
s  w 


(1.18) 


(1.19) 


(1.20) 


(1.21) 


Consider  first  the  homogeneous  form  of  Equation  (1.16),  which  is 


(D4  +  bD3  +  cD2  +  dD  +  f)x  =  0  (1.22) 

If  we  assume  a  solution  of  the  form 

xH(t)  =  Ceat  (1.23) 

then  substitution  of  Equation  (1.23)  into  Equation  (1.22)  yields 

(a4  +  bcx3  +  ca2  +  dd  +  f)Ceot  =  0  (1.24) 

Assuming  that  Ceat  £  0,  Equation  (1.24)  yields 

a4  +  ba3  +  Ca2  +  da  +  f  =  0  (1.25) 

The  solution  of  Equation  (1.25)  is  obtained  in  Appendix  J. 

Having  solved  Equation  (1.25),  we  write  Equation  (1.23)  in  the  form 

4  o.t 

x  (t)  =  T  C  e  J  (1.26) 

H  j=l  J 

We  now  return  to  Equation  (1.16),  which  can  be  written  in  factored  form 

[(D-a1)(D-a2)(D-a3)(D-a4)]x  =  gDp  (1.27) 


The  particular  solution  to  Equation  (1.27)  is  developed  in  Appendix  K. 

The  complete  solution  to  Equation  (1.16)  is  the  sum  of  the  homogeneous 
solution,  given  by  Equation  (1.26)  and  the  particular  solution,  given  by 


-% 
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(1-28) 

(1.29) 

(1.30) 

(1.31) 

(1.32) 

(1.33) 

(1.34) 
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t 


€ 


C> 


p(0)  +  OjP(0)  +  OjP(0) 


o 


(1.35) 


When  p(t)  is  a  step  pulse,  i.e.,  when 

p(t)  =0  (t<0) 

=  P  (t>0) 


(1.36) 


then 


p( o)  =  p(0)  =  p(0)  =  p(0)  =  0 


(1.37) 


so  that  for  a  system  initially  at  rest.  Equations  (1.32),  (1.33),  (1.34) 
and  (I .35)  reduce  to 


1 

1 

1 

1 

(  0  N 

“1 

a2 

“3 

a4 

c2 

0 

2 

2 

2 

2 

- 

> 

al 

a2 

a3 

a4 

C3 

0 

3 

3 

3 

3 

“1 

a2 

°3 

a4 

^  C4  y 

(1.38) 


The  coefficient  matrix  in  Equation  (1.38)  is  nonsingular,  provided 

*  «j  (1«)  (1-39) 

because  its  determinant  is  a  Vandermonde  determinant  [Bellman  (1960:186)], 
and  therefore  it  must  be  that 

(1.40) 

In  this  case  no  homogeneous  solution  is  required  to  satisfy  the  initial 
conditions,  and  therefore  the  particular  solution  is  the  complete 
solution.  Equation  (1.28)  yields 


ci  -  c2  -  c3  =  C4  -  0 


A  i  ait 
M1)  =  xps(t)  =  gP  A.  ~  e  3 


j=i  ’j 


(t>o) 


(1.41) 
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so  that 


*s(t)  =  gP  L  ir 
s  j-i  'j 


t  ^  *°jt 


(t>0) 


(1.42) 


4  a.  a -t 

(t)  =  gp  Z  ir  e  J 

j=i  "j 


(t>0) 


(1.43) 


From  Equations  (1.41),  (1.42)  and  (1.43)  we  obtain 


M°+)  ■=  I  ~ 


j-i  ’J 


lM“2-a3Ma2~a4; 

1  „  1 
(a3-al )  (a3-a2^  (a3-a4)  ( a4_al )  ( a4~a2  )  ( 0,4~“3  ) 

2  2 

l  [«2  “  (a3+a4^°2  +  a3a4-l  “  t°i  "  (°3+a4^ai  +  °3a4^ 


a^-a2 )  (al“a3^ (ai-cV (a2“  “3^ (“2"  “4 

2  2 
l  [°4  -  («i+“2^a4  +  ala2^  “  ^-“3  “  (al+a2^a3  +  “la2^ 


a3-a4 )  ( a3~ai ) ( a3_a2  > ( a4~ai ) ( a4~a2 


(a3+a4) (al_a2^  “  (al+a2^al~c,2^ 

[al-a2 )  (al-a3 )  (ai-a4 )( a2— a  3 )  (a2-ct4J 


(al+a2 ^ ^a3-a4^  “  ^a3+a4^a3~a4^ 


'a3-a4'  'a3-otl '  ^a4-al'  'a3-ot2  'a4_a2; 


(1.44) 
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Cf 


( “3-°>4 )  ( 03_“i )  ( a3-a2  > '  a4-ai )  ( °4'a2 ' 

(al+a2^ ^al_a2 ^a3a4  ~  (a3+ot4 ) ( “i- a2 ^ al a2 

( a i~ a 2  ^  ( 01 1-01 3 )  ( a i_ a 4 )  ( a2— a3 )  (a2“a4) 

^a3+a4^a3_a4^ala2  _  (0tl+a2^a3~0t4^<:,3a4 
(a3_a4^  (a3-“i)  (a3-a2^a4"al^a4"a2^ 


=  0  (1.46) 

A  numerical  example  shows  that  two  roots  of  Equation  (1.25)  are  real  and 


the  remaining  two  are  a  complex  conjugate  pair:  All  real  quantities  are 

negative. 

=  -m 

(1.47) 

a2  =  -n 

(1.48) 

a3  =  -q  +  is 

(1.49) 

a4  =  -q  -  is 

(1.50) 

Equation  (1.41)  can 

therefore  be  written  in  the  form 

x$(t)  =  gP 

1  P-mt  +  1  P~nt  +  1  P-qt+ist  +  1  p-qt-ist 

L  *1  "2  *3  *4  J 

(1.51) 

Now 

i 

e  =  cos  st  +  i  sin  st 

(1.5 2) 

e-lst  =  cos  st  -  i  sin  st 

(1.53) 

so  that  Equation  (1.51)  can  be  written  in  the  form 


*• 
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x.(t)  =  gP  e_mt  +  -i  e-nt 

S  *2 


*  (-UJA 

V’3  "4  / 


e  qt  cos  st 


+  ,(4'^)e‘qt  sinst]  (I-54) 


x  (t)  =  gP(Ae-mt  +  Be“nt  -  D1e_qt  cos  st  +  D2e~qt  sin  st)  (1.55) 


where 


A__l=  1 _ 1 

-  v1  ~  (n-m)[(q-m)2  +  sZ] 


(1.56) 


B  =  —  = 


*2  -  (a2-a1)(a2-a3)(o2-a4)  “  (n_m)[(q_n)2  +  $2] 


(1.57) 


D.  =  -  —  +  —  [See  last  two  terms  of  Equation  (1.44).] 

l  ir3  w4 

(«1+«o)  -  (<*-}+°a) 


l  a3-ai > \ a4-al t V a3-a2 ' 1 a4_a2 ; 


_ m  +  n  -  2q 

[<nv-q)z  +  s^Jttn-q)2  ♦  s‘ 


(1.58) 


D  _  i  _1  _  JL  ,  i  i 

2  ir3  ir4  (<*3-  ap(a3-a2)(a3-a4J 


l  a4-ai / 1 a4~a2 / v a4~a3  J 


2  2 
i  [a4  -  (apa2 )a4  +  aia2^  +  ^a3  “  ^al+a2^a3  +  ala2^ 


a3"°4 


l0l3-al  ^°3_a2' 'a4-al'  'a4~a2' 


i  C ( c,3+c*4 )  "  (al+a2 ^a3+a4^  +  2ala2^ 
P*4  (a3-al  )  (a4_al )  (°3“a2H®4-“2^ 


i  [2q2  -  2s2  -  (m+n) ( 2q)  +  2mn' 


2lS  [(m-q)2  +  s2][(n-q)Z  +  s2] 
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2  2 

_  mn  -  (m+n) q  +  q  -  s 
'  s[(m-q)*  +  s2][(n-q)2  *  s2] 

Equation  (1.55}  can  be  simplified  by  a  trigonometric  substitution. 


If  we  set 

R  =  Vd  \  +  0* 

-1  °2 
3  =  tan  fr 

U1 

then  Equation  (1.55)  can  be  written  in  the  form 

x$(t)  =  gP[Ae_mt  +  Be-nt  -  Re'qt  cos(st+s)] 

For  the  case  at  hand.  Appendix  L  gives 

m  =  0.115726598  X  10-4  MSEC-1 
n  =  0.538700000  X  107  MSEC-1 
q  =  0.462842340  X  10-6  MSEC-1 
s  =  0.425932181  X  10°  MSEC-1 


(1.59/ 


(1.60) 

(1.61) 

(1.62) 
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(ft 


'ft 


'ft 


'ft 


T 


g  =  1.508  X  10° - — - * 

LB(MSEC)"3 

so  that 

(m-q)2  +  s2  =  0.1814182229  MSEC"2 
(n-q)2  +  s2  =  2.9019769  X  1013  MSEC"2 

and  therefore 

A  =  - ^—5 - 5-  =  1.023227294  X  10~6  (MSEC)3 

(n-m)[(q-mr  +  s2] 

B  =  - L_ -  =  -6.396745516  X  10“21  (MSEC)3 

(nv-n)[(q-n)2  +  s2] 

D,  = -  m  *  n  -  2q  -  ,  1.023227293  X  10"6  (MSEC)3 

1  [(m-q)2  +  s2][(n-q)2  +  s2] 

D0  = - mn  -  (^n)q  +  <£_-_s2  =  2.660848297  X  10_11(MSEC)3 

2  s[(m-q)2  +  s2][(n-q)2  +  s2] 

R  =  Vd2  +  D2  =  1.023227293  X  10~6  (MSEC)3 

-1  ^2  -3 

B  =  tan  1  yf  =  1.489942369  X  10  J  DEG 

U1 

so  that  Equation  (1.62)  gives 

^11  =  )  543026759  X  10'6  e'0’ 115726598  X  *** 

-  9.646292238  X  10“21  e-°-538700000  X  ^ 

-  1.543026758  X  KT6  e-°-462842340  X  10'6‘ 

*  [cos(0.425932181t  +  2.600440001  X  10"5)]  FT/LB 


'ft 


L. 
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Differentiation  of  Equation  (1.62)  yields 

x s ( t)  =  gP[-mAe-mt  -  nBe~nt  +  qRe-qt  cos(st+e) 

+  sRe-c,t  sin(st+B)]  (1.63) 

Equation  (1.63)  can  be  simplified  by  another  trigonometric  substitution. 


If  we  set 

r  =  V  q2  +  s2  =  0.425932181  MSEC-1  (1.64) 

t  =  tan  -1  q/s  =  6.225759401  X  10-5  DEG  (1.65) 

then  Equation  (1.63)  can  be  written  in  the  form 

xs(t)  =  gP[-mAe-mt  -  nBe-nt  +  rRe-qt  sin(st+e+0)]  (1.66) 

where 

S  +  =  1.552199963  X  10-3  DEG 

so  that  Equation  (1.66)  gives 

=  -1.785692374  X  10-11  e"0-115726598  X  10  t 
+  5.196457629  X  10-14  g"0-538700000  x  1()7t 
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+  6.572247524  X  10“7  e"0*462842340  X  10  6t 

*  [sin(0.425932181t  +  2.709100001  X  10"5)]  (1.67) 

FT/ (LB  MSEC) 

Differentiation  of  Equation  (1.66)  yields 

xs(t)  =  gP[m2Ae  +  n2Be  -  qrRe  sin(st+e+0) 

+  srRe“qt  cos(st+&+0)] 

-  gP[m^Ae-mi'  +  n2Be“nt  +r2Re-qt  cos(st+e+20)]  (1.68) 

where 

e  +  20  =  1.614457557  X  10-3  DEG 
so  that  Equation  (1.68)  gives 

=  2.066521035  X  10"16  e"°' 115 726598  X  10  4t 

-  2.799331725  X  10-7  e-°* 538 700000  x  1C)7t 
+  2.799331722  X  10-7  e~° • 462842340  x  10 

*  [cos(0.425932181t  +  2.817760001  X  10"5)]  (1.69) 

FT/LB(MSEC)2 

The  impulse  response  function  can  be  obtained  from  Equation  (1.62),  or 
from  the  fact  that  an  impulse  is  the  successive  application  of  two  step 
pulses. 
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<# 


at 


i°(t) 


•• 


h(t) 


xs(t)  -  X$(t-At)  ^  [x$(t)]At 
PAt  =  PAt 


_  *s(t)  =  g[-mAe_mt  -  nBe 

'  P 

=  -  1.785692374  X  10"11 
+  5.196457629  X  10-14  e 


nt  +  rRe  sin(st+e+0)] 

-0.115726598  X  10~4t 
e 

-0.538700000  X  107t 


(1.70) 


♦  6.572247524  X  KT7  e-°'462842340  *  10~6t 

•  [sin(0.425932181t  +  2.709100001  X  10~5)]  (1.71) 

FT/ (LB  MSEC) 

The  model  associated  with  Equations  (1.5)  or  (1.6)  is  identical  to 
that  associated  with  Equations  (H.5)  or  (H.6),  shown  following  Equation 
(H.27),  except  that  the  pistons  have  mass. 
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I 


A  suddenly  applied  force  p(0+)  will  produce  no  initial  acceleration  of 

mass  because  there  will  be  no  initial  relative  velocity 

y(0+)  -  x(0+),  and  therefore  no  initial  dashpot  force  to  cause  such  an 

acceleration.  For  awhile  the  mass  M  will  drag  the  mass  M  along  with 

it,  but  when  the  system  comes  to  rest  the  mass  will  return  to  its 

initial  position  and  the  spring  K  will  carry  the  entire  final  load. 

w 
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APPENDIX  J 


Solution  of  a  Fourth  Order  Algebraic  Equation 

Consider  the  fourth  order  algebraic  equation  with  constant  coefficients 
a4  +  ba3  +Ca2  +  da  +  f  =  0  (1.25) 

The  solution  of  Equation  (1.25)  is  classic,  but  although  involved,  the 
resulting  equations  need  not  be  as  useless  as  is  commonly  believed 
[Sokolnikoff  and  Sokolnikoff  (1941:91)]. 

To  eliminate  the  cubic  term,  let 


°  -  2  "  4 


(J.l) 


so  that 


a  = 


,4  „3,  .  3  2,2  1  .3  .  1  . 

z  -  zb  +  ^zb  -  jg  z  b 


(J.2) 


3  3  3  2,  .  3  .2  1  .3 

a  .  z  -?zb+Hzb  -^b 


(J.3) 


2  2  1  .  ,  1  ,2 
a  =  z  -  ?  zb  +  T6  b 


(J.4) 


Substituting  Equations  (J.l),  (J-2),  (J.3)  and  (J.4)  into  Equation  (1.25) 
yields 

,4  3.  .  3  2,2  1  ,3  ,  1  ,4 

z  -zb+gzb  -jg2b  256  b 


,  ,/  3  3  2,  .3.2  1  ,3. 

+  b(z  ~  7T  z  b  +  zb  -nb) 


64 


+  c(z2  -  i  zb  +  b2) 


+  d(z  -  i  b)  +  f 
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i 


ft 


% 


=  z4  +  (-  |  b2  +  c)z2  +  b3  -  \  be  +  d)z 
+  256  b  +  T6bc-4bd  +  f)=0 


(0.5) 


or 


rz  +  s 


0 


(0.6) 


where 


(0.7) 


13  1 

r=_gb  ~IbC  +  d 


(0.8) 


s 


b2c 


-  i  bd  +  f 


(0.9) 


Equation  (J.6)  can  be  expressed  as  the  product  of  two  second  order 
polynomials,  having  their  linear  terms  equal  in  magnitude  and  opposite  in 
sign,  of  the  form 

[(z-1)2  -  (m+n)2][(z+l)2  -  (m-n)2] 

=  (z-l)2(z+l)2  -  (z-l)2(m-n)2  -  (z+l)2(frd-n)2 
+  (m+n)2(n>-n)2 

=  (z2-!2)2  -  (z2-21z+12) (m2-2mn+n2) 

-  (z2+21z+l2)(m2+2mn+n2)  +  (m2-n2)2 
=  (z4-212z2+l4)  -  (m2-2mn+n2)z2 

+  (21m2-41mn+21n2)z  -  (l2m2  -  212mn+l2n2) 

-  (m2+2mn+n2)z2  -  (21m2+41mn+21n2)z 

-  ( l2m2+212mn+l2n2)  +  (m4-2m2n2+n4) 

=  z4  -  2( 12+m2+n2)z2  -  8  lmn  z 

+  [(l4+m4+n4)  -  2(  l^m^+m^n^+n^^)]  _  o 

(J.10) 
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Comparison  of  Equations  ( J . 6 )  and  (J.10)  shows  that 


-2(12  +  m2  +  n2)  =  q 

(0.11) 

-8  Inn  =  r 

(0.12) 

(1  +  m  +  n  )  -  2(1  m  +  m  n  +  n  1  ]  =  s 

(0.13) 

so  that 

l2  +  m2  +  n2  =  - 

(0.14) 

lV  ♦  n,2n2  -  n2!2  .  '  s  .  0^*1 

(0.15) 

,222  r2 

1  m  "  =  64 

(J.16) 

Equations  (J.14),  (J.15)  and  (J.16)  define  the  coefficients 

of  a  cubic 

2  2  2 

equation,  of  which  1  ,  m  and  n  are  the  roots. 

(0.17) 

When  1,  m  and  n  have  been  found  by  solving  Equation  (J.17), 

being  careful 

to  ensure  that 

Imn  =  - 

( J.18) 

then  the  initial  form  of  Equation  (J.10)  can  be  written  in 

the  form 

(z-l-m-n)  (z-l+nr+n)  (z+l-m+n)  (z+l+m-n)  =  0 

(0.19) 

and  therefore  the  four  roots  are 

Zj  =  1  +  m  +  n 

( J.20) 

c 

1 

e 

i 

a 

CM 

N 

(0.21) 

z-j  =  -1  +  m  -  n 

(J.22) 

z4  =  -1  -  m  +  n 

(J.23) 

Equation  (J.17)  can  be  written  in  the  form 
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(0.24) 


k3  +  Bk2  +  Yk2  +  «  =  0 


where 


y  = 


(0.25) 

q2  -  4s 

16 

(0.26) 

2 

u|S 

1 

II 

(0.27) 

In  order  to  eliminate  the  second  order  term  in  Equation  (J.24),  we  set 


k  =  w  - 


(J.28) 


so  that 


?  3  2  12  13 

k  =  w  -  ew^  +  “  27  B 


k2  =  w2  -  |  BW  +  |  B2 


(0.29) 

(0.30) 


Substituting  Equations  (J.28),  (0.29)  and  (0.30)  into  Equation  (0.24) 
yields 

(w3  -  bw2  +  ^  b2w  -  yj  B3)  +  b(w2  -  4  0w  +  ^  b2)  +  y(w  -  4)  +  6  =  0 


or 


W3  -  (4  -  Y)W  -  (*J  -  ^  -  6)  =  0 


(J.31) 


or 


w°  -  Pw  -  Q  =  0 


(J.32) 


where 
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To  obtain  a  solution  of  Equation  (J.32),  we  assume  that 


w  =  A  +  B  (J.35) 


so  that 


w3  =  A3  +  3A2B  +  3AB2  +  B3 
=  A3  +  B3  +  3AB(A  +  B) 

=  A3  +  B3  +  3ABw 

(J.36) 

Equation 

(J.36)  can  Lt  written  in  the  form 

w3  -  3ABw  -  (A3  +  B3)  =  0 

(J.37) 

Comparison  of  Equations  (J.32)  and  (J.37)  shows  that 

3AB  =  P 

(J.38) 

A3  +  B3  =  Q 

(J.39) 

Equation 

(J.38)  yields 

A3B3  =  (£)3 

(J.40) 

Thus  the 

3  3 

sum  and  product  of  A  and  B  are  specified,  which  means 

that 

they  are 

the  roots  of  the  quadratic  equation 

Z2  _  (a3  +  B3)£  +  (A3B3)  =  0 

or 

52  -  Q5  +  (y)3  =  0 

(0.41) 

Thus 

El  -  A3  .  §  ♦  V(f)2  -  <^)3 

(0.42) 
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3 


s2  =  B3  =  §  -  V<§)2  -  (j)3  (J.4 

If  we  set 

R  =  (Q)2  -  (^)3  (J.44) 


then  Equations  (J.42)  and  (J.43)  can  be  written  in  the  form 


A3  .  | 

(J.45) 

B3  =  |  ^ 

(J.46) 

When  R<0,  i.e. ,  when 

<£)3 » <$>* 

(J.47) 

then  Equations  (J.45)  and  (J.46)  take  the  form 

A3  .  §  ♦  i  V(|)3  -  (§)2 

(J.48) 

B3=§-  i  V(|)3  -  (f)2 

(J.49) 

3  3 

In  this  case,  A  and  B  can  be  represented  in  the  complex  plane  as 
shown  below. 


Thus  we  can  write 


a3  ,  Vtf)5  e13e 

( J .  50 ) 

B3  =  V^)3  e~'38 

(J.51) 

where 

-i  1 

3e  =  cos  - - - 

V  <^)3 

(J.52) 

so  that 

Aj  -Vf  e**  (j  =  1,  2,  3) 

(J.53) 

Bk  »  V|  e  (k  =  1,  2,  3) 

(J.54) 

where 

i  -1  I 

0,  -  e  -  -  cos  - 

V<>)3 

(0.55) 

*2  ■ 8  *  H 

(0.56) 

i  2w 

^3  =  e  J 

(0.57) 

Now  Equation  (J.38)  requires  that 

AB  =  j 

(J.58) 

so  that  the  choice  of  Aj  and  in  the  relation 

wj  -  Aj +  Bk 

( J .59) 
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► 

> 

f 

► 

t 

must  be  such  that  and  are  complex  conjugates. 

Therefore 

it  must 

> 

be  that  k  =  j,  so  that 

!• 

IT  i^i 

w j  =  V  ^  (e  J  +  e  J) 

• 

=  2  V|"  cos  t>.  (j  =  1,  2, 

3) 

(J.60) 

The  above  three  roots  are  all  real. 

When  R>0,  i.e.,  when 

• 

<|)2  >  (|)3 

( J.61) 

3  3 

then  both  A  and  B  ,  as  defined  by  Equations  (J.45) 

and  ( J.46) 

are 

• 

real. 

If  we  set 

• 

°i =  1 +  V* 

(J.62) 

°2  =  §  "  V** 

( J . 63 ) 

• 

and 

Ao  =  !oxl1/3  sgn  ax 

( J .64) 

<• 

Bq  =  ia2|1/3  sgn  a2 

(J.65) 

then 

• 

Aj  =  Aoe  ° 

( J.66) 

i*k 

Bk  =  Boe 

( J.67) 

• 

where 
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01 

.  ■  ^ 

\  =  o 

(J.68) 

* 

±  2ir 

h  =  ~J 

(J.69) 

A  2* 

^3  =  J 

(J.70) 

• 

Now 

Vo  -  V  ^  ^  *  <5>3  ■  T 

( J .  7 1 ) 

• 

which  means  that  the  choices  of  Aj  and  B ^  in  the  relation 

wj  =  Aj  +  Bk 

(J.59) 

ft 

must  be  such  that  and  are,  again,  complex  conjugates,  and 

therefore 

W1  -  A1  *  B1  -  Ao  +  Bo 

(J.72) 

• 

.  2ir  ■  2ir 

1  “1-T 

w9  «  A,  +  B,  =  A  e  +  B  e 

2  2  3  o  o 

• 

-  -  7  <Ao  +  Bo>  *  1  *  <\  -  Bo> 

(J.73) 

_  i  i  % 

w9  =  A  +  B9  =  A  e  +  B  e 

3  3  2  o  o 

'ft 

=  "  7  (Ao  +  Bo)  ‘  1  "T  (Ao  ‘  V 

(J.74) 

In  summary,  to  solve  Equation  (1.25)  we  do  the  following: 

ft 

t  Cv<Ms*Mw> 

'  MsMw 

(1.17) 

'  ft 

.  Vw+Vs 

C  '  Vw 

(1.18) 

f  ft 

I. - 
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(1.19) 


d  = 


Cv(Ks  +  Kw> 
MsMw 


f  = 


KsKw 

MsMw 


then 


then 


then 


Cv 

9  '  MsM„ 


3  .2  .  r 

q  -  -  8  b  +  c 


r  @  1  b3  _  1  bc  +  d 


s  =  "  235  b4  +  T3  fa2c  ‘  \  bd  +  f 


“I 


q  -  4s 
Y  =  1’6 


5  -  "64 


P  r 


0  _  By  _  2e 
Q  =  1  ^7 


-  5 


then 


(1.20) 

(1.21) 

(0.7) 

(0.8) 

(0.9) 

( J .25) 

(0.26) 

(0.27) 

(0.33) 

(0.34) 
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r  =  qr  -  (§)a 


(J.44) 


If  R<0: 


0  =  "j  cos 


(J.55) 


2 cos  o  - 


=  2V4  cos(o  +  — 5")  - 


(0.75) 


(0.76) 


=  2V?-  cos(e  -  - 


(0.77) 


kj  =  MAX(k1,  k2,  k3) 


(0.78) 


ep  =  cos 


-1  /3kI  +  e 


(0.79) 


kII  =  2  cos(ep  +  %  "  1 


(0.80) 


kl II  “  2  ^3"  cos(ep  I 


(0.81) 


n  =  ^  cos  1(sgn  kn) 


(0.82) 


1  =  - 


r  cos  2« 


8I  kI 1 1  I  kII 1 1 


(0.83) 


=  | kjj j  Z  (cos  ft  +  i  sin  ft) 


(0.84) 


=  0  or  it  /  2 
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1 


ft 


n  =  (cos  ft  +  i  sin  ft) 


( J .85) 


The  rationale  for  Equations  (J.83),  (J.84)  and  (J.85)  is  explained 
following  Equation  (J.94). 

If  R>0: 


°1  =  7  +  V* 


°2  '  2 


-V* 


A0  =  I  1 1/3  sgn  oj 


80  =  U2I1/3  sgn  a2 


»  -  j&Ao  +  Bo>  +  fJ2  +  l(Ao  -  Bo>2 
7(Ao  +  V  +  T 


1/4 


1  -1 

ft  =  cos  ^ 


{_ 


1  =  - 


8P 


( J.62) 
(J.63) 
(J.64) 
(J.65) 
( J.86) 

( J.87) 

(J.88) 


m  =  p(cos  ft  +  i  sin  ft)  (J.89) 

n  =  p ( cos  SI  -  i  sin  ft)  (J.90) 

The  rationale  for  Equation  (J.88)  is  explained  following  Equation  (J.94). 

Equations  (J.72),  (J.73),  (J.74),  (J.86)  and  (J.87)  can  be  represented 
graphically  as  shown  below. 


>  i 
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t _ 


rf  is  the  complex  conjugate  of  mu  (mirror  image  across  the  x-axis). 
Finally,  we  have 


ojsl  +  m+  n-  ^- 


(J.91) 


a  2  =  1  -  m  -  n  -  ■y 


(J.92) 


a3  =  +  ro  -  n  -  ^ 


(J.93) 


a4  =  -l-m+n-^ 


(J.94) 


Equations  (J.83)  and  (J.88)  ensure  that  1  will  be  real,  and  that  Equation 
(J.18)  will  be  satisfied.  The  correctness  of  the  assumption  that  1  is 
real  can  be  verified  as  follows: 

Notice  that  Equation  (J.16)  requires  that  the  product  of  the  three 

2  2  2 

roots  of  Equation  (J.17),  i.e.  1  m  n  ,  be  a  nonnegative, 
real  number. 
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?  ?  2 

If  the  three  roots,  1,  m  and  n  ,  are  all  real,  then  either  one 

is  nonnegative  and  the  other  two  are  negative,  or  all  three  are 

2  2 

nonnegative.  In  either  case,  m  and  n  have  the  same  sign, 
which  is  assured  by  Equations  (0.84)  and  (J.85). 

If  one  root  is  real  and  the  other  two  are  complex,  then  the  complex 
roots  must  be  complex  conjugates.  Since  the  product  of  the 
complex  conjugate  roots  is  nonnegative,  the  real  root  must  also  be 
nonnegative. 

Thus,  we  can  always  assume  that 

l2  >  0  (0.95) 

and  therefore  that  1  is  real. 

The  coefficients  (invariants)  of  a  fourth  order  polynomial  equation 
will  be  real  if  and  only  if  complex  roots  occur  only  as  complex  conjugate 
pairs. 


*1  =  al  +  °2  +  a3  +  a4  =  ^al+a2^  +  ^a3+a4^ 

(J.96) 

l2  =  ala2  +  ala3  +  ala4  +  a2a3  +  a2“4  +  a3a4 

=  0^2  +  (ai+a2 )  (a3+a4 )  +  a3a4 

(0.97) 

T3  =  a2a3a4  +  ala3a4  +  aia2°4  +  ala2°3 

=  “304(0^02)  +  ala2^0l3+a4^ 

(J.98) 

I4  =  °ia2a3a4  =  (aia2^  ^a3a4^ 

(0.99) 
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APPENDIX  K 


Solution  of  a  Fourth  Order  Linear  Ordinary  Differential  Equation  With 
Constant  Coefficients 


Consider  the  fourth  order,  linear  ordinary  differential  equation  with 


constant  coefficients 

(D4  +  bD3  +  cD2  +  dD  +  f)x  =  gDp  ( 

which  can  be  written  in  the  factored  form 

[(D-a1)(D-a2)(D-a3)(D=a4)]x  =  gDp  ( 

where  «2,  a3,  and  a4  are  the  roots  of  the  equation 

a4  +  ba3  +  Cci2  +  da  +  f  =  0  ( 

The  particular  solution  to  Equation  (1.27)  can  be  written  in  the  form 


(1.16) 


(1.27) 


(1.25) 


Xp(t}  =  (D-ai)(D-a2)(D-a3)(D-a4)  gDp 


(K.l) 


Now  if  we  set 


_A_  +  B 


3-a^ ) ( D-a2)  D-a^  D-a2 


(K.2) 


then  it  must  be  that 


and  therefore  that 


AD  -  Aa2  +  BD  -  Boj  =  1 


(K.3) 


B  =  -A 


(K.4) 


so  that 


and  therefore 


A(a^  -  a2)  =  1 


(K.5) 


al  ”  a2 


( K .  6 ) 
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2  ~  al 


(K.7) 


Thus,  we  have 


'“l  2' '  al'  ^  2  1M  2J 


(K.8) 


If  we  set 


•i  -  m  -  Ru  <''«> 


then  for  distinct  roots. 


D  -  a-  =  D- 

l  i 


_1 _ 1  X  1 

°iDrVi  Vj 


(K.9) 

(K.10) 


(K.ll) 


We  then  have 


)1D2D3  ^R12D1  R21D2^^ 


+  _X_)  ♦  1(1  +  _1_) 

R12  R13°l  R31D3  R21  R23D2  R32D3 


_ 1 _ +  _ 1 _ +  J_(  J _ l_)_i 

R12R13D1  R21R23D2  R12  R31  R32  °3 


1  +  1  +  R32  ~  R31 

R12R13D1  R21R23D2  R12R31R32D3 


_ 1 _  +  _ L__  +  _ L_ 

R12R13D1  R21R23D2  R31R32°3 


(K.12) 


and,  in  general. 


Dl°2°3  * 


_ l  -L 

■  dn  j-1  'jDj 


(K.13) 
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where 


11  j  =  RjlRj2Rj3  Rj,j-lRj,j+l  '•*  RjN 
Now  the  particular  solution  for  is,  from  Equation  (H.15), 

J 


(K. 14) 


i  t  .  a.(t-x) 

IT  £*&-*'$*  *  dl  (t>0> 

3  ^ 


(K. 15) 


Equation  (K.13)  therefore  yields 


xp(t)  '  d1D2D3D4  ^  "dt^  ' 


4  it.  a.(t-x) 

9  I  — J  77  e 
j=l  ’j  0  dX 


(t>0) 


(K.16) 
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APPENDIX  L 

Numerical  Solution  of  the  Characteristic  Equation 

Equation  (1.25)  will  be  solved  for  one  cubic  foot  of  dry  sand, 
constrained  to  undergo  one  dimensional  compression  and  loaded  by  a  100  PSI 
airblast. 

From  [Terzaghi  and  Peck  (1967:28)]  we  obtain 
Ws  =  115  LBS 


M  -  _  3  5714  SEC 

"s  ”  ~  fT^ 


n  =  0.30 


VA  =  Vv  =  0.30  FT3 


and  since  [McGuire  (1968:181)] 


Ya  =  0.07651 
rt  FT-3 


we  have 


WA  =  (0.07651) (0.30)  =  0.022953  LBS 


u  0.022953  ,  V  ,„-4  LB  SEC‘ 

"w  -  32.2  -  -  7-1Z8  X  10  FT- 


From  [Bowles  (1977:269)]  we  have 


Ks  =  600,000  m 


and  if  we  assume  isothermal  compression,  we  obtain 
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d£  _  dV 
P  ‘  “  V 


(L.l) 


so  that 


dp  =  -  ^  dV 


(L.2) 


A  dp  =  £y-*A  dy 


(L.3) 


Kw  .  ja!  ■  .  4M00  LBS 


Assuming  an  effective  permeability,  as  given  by  Equation  (4.3),  of 


B1  “  5,054  X  10  n  SEC 


(0.3048) 


=  2.605  X  10 


-7  FTh 
LB  SEC 


we  have 


Cy  =  ~(1  FT3)  =  3,839,249  LBF3E- 


Thus,  we  obtain 


CV(MS  +  Mw) 


(3.839  X  10°) (3.5714  +  0.0007] 
(3.571)  0.0007128) 


=  5.387  X  109  SEC-1 


MsKw  *  MwKs  _  (3.571) (48,000)  +  (0.0007128) (600,000) 


[0.0007128] 


=  6.751  X  107  SEC"2 
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VKs  *  Kw>  (3.839  X  106)  (600,000  ♦  48,000 ] 


[3. 571 ) (0. 0007128) 


=  9.773  X  1014  SEC"3 


*  KsKw  (600,000) (48,000)  ,  01  v  ml3  ccr-4 

f  =  o~  =  tTJTTjfotmrim =  SE 

j  n 


3.839  X  10c 


Changing  from  SEC  to  MSEC,  we  obtain 
b  =  5.387  X  106  MSEC-1 
c  =  6.751  X  101  MSEC"2 
d  =  9.773  X  105  MSEC"3 
f  =  1.131  X  101  MSEC'4 


=  1.508  X  10“ 


LB  SEC' 


g  =  1.508  X  10L 


LB  MSEC' 


The  following  double  precision  FORTRAN  computer  program  was  used  to 
calculate  the  roots  of  Equation  1.25.  The  only  differences  between  the 
program  instructions  and  the  equations  given  in  Appendix  J  are  due  to  lack 
of  double  precision  functions  for  the  inverse  sine  and  cosine. 

PROGRAM  AIR(INPUT, OUTPUT) 

DOUBLE  PRECISION  B,C,D,F,Q,P,S, B1 ,C1 ,D1 ,PP ,QQ,RR ,TH 

1  C0S1 ,C0S2 ,C0S3,THP ,Y, SQL , SQM,SQN, RM,RN,XL, SM,YM, 

2  XN,YN,S1,S2,EX,A0,B0,RH0, 

3  OMEGA, XR1,YR1,XR2,YR2,XR3,YR3,XR4,YR4 
B  =  0.5387D+07 

C  =  0.6751D+02 
D  =  0.9773D+06 
F  =  0.1 1310+02 
P I =4 . 0*DAT AN ( 1 . 0D0 ) 

Q=-0.375*B**2+C 

R=0 . 125*B**3-0 . 5*B*C+D 

S=-(3.0*B**4)/256.+(C*B**2)/l6.0-0.25*B*D+F 
B 1 =Q / 2 . 0 
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C1=(Q**2-4.0*S)/16.0 

Dl=-R**2/64. 

PP=(Bl**2)/3.0-Cl 

QQ=( Bl*Cl)/3.0-(2.0*Bl**3)/27.0-Dl 
RR=  (QQ/2.0)**2-(PP/3.0)**3 
IF(RR.GT. 0.0)60  TO  10 
TH=QQ / ( 2 . 0*DSQRT (-RR ) ) 

TH=DATAN(TH) 

TH= (P I / 2 . 0— TH ) / 3 . 0 
C0S1=DC0S(TH) 

C0S2=DC0S(TH+(2.*PI)/3.0) 

C0S3=00CS(TH-2.0*P I )/ 3.0) 

C0S1=DMAX1(C0S1 ,C0S2 ,C0S3) 

T  H  P=D  SQRT  ( 1 . 0-CO  SI*  *2 )  /  CO  SI 
THP=DATAN(THP) 

C0S2=0C0S(THP+(2.0*PI )/3.0) 

C0S3=D0C S ( THP-( 2 . 0*P I ) / 3 . 0 ) 

Y=2.0*DSQRT (PP/3.0) 

SQL=Y*C0Sl-Bl/3.0 
SQM=Y*C0S2-Bl/3.0 
SQN=Y*C0S3-Bl/3.0 
RP4=DSQRT ( DABS(  SQM) ) 

RN=DSQRT ( DABS( SQN ) ) 

0MEGA=0 . 5*( P I / 2 . 0-2 . 0*DATAN ( DS I GN ( 1 . ODO , SQM+SQN ) ) ) 
XL=DSQRT(SQL) 

XM=RM*DCOS( OMEGA) 

YM=RM*DSIN(OMEGA) 

XN=RN*DCOS ( OMEGA ) 

YN=RN*DSIN (OMEGA) 

GO  TO  20 

10  S1=QQ/ 2 . O+DSQRT ( RR ) 

S2=QQ/2.0-DSQRT(RR) 

EX=l./3. 

A0= ( ( DABS ( SI ) ) **E  X ) *DSI GN ( 1 . ODO , SI ) 

B0=( ( DABS( S2 ) ) **EX ) *D  SI  GN ( 1 . ODO , S2 ) 
RH0=((0.5*(A0+B0)+Bl/3.0)**2+0.75*(A0-B0)**2)**.25 
0MEGA=( A0+B0+2.0*Bl/3.0)  /  (DSQRT (3.0D0)*(A0-B0) ) 
OMEGA=DATAN( OMEGA) 

OMEGA=O.5*(PI/2.O+0MEGA) 

XL=-R/(8.0*RH0**2) 

XM=RH0*DC0S( OMEGA) 

YM=RHO*DSIN (OMEGA) 

XN=XM 

YN=-YM 

20  XRl=XL+XM+XN-B/4.0 
YR1=YM*YN 

XR2=XL-XM-XN-B/4.0 

YR2=-YM-YN 

XR3=-XL+XM-XN-B/4.0 

YR3=YM-YN 

XR4=-XL-XM+XN-B/4.0 

YR4=-YM+YN 
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100 


PRINT  100,B,C,D,F,XR1 ,YR1,XR2,YR2,XR3,YR3,XR4,YR4 
FORMAT(10X,2D16.9) 

STOP 
END 

The  roots  are 

oj  =  -0.115726598  X  10-4  MSEC'1 
o2  =  -0.538700000  X  107  MSEC-1 
a3  =  -0.462842340  X  10-6  +  0.42593218H  MSEC-1 
a4  =  -0.462842340  X  1 O-6  -  0. 425932181 i  MSEC-1 
The  invariants  are 

*1  =  “l  +  a2  +  a3  +  “4  =  “5,387,000 

12  =  0^2  +  (aj+a.,)  (t*3+a4)  +  a3a4  =  67.51 

13  -  (a1+a2)a3o4  +  «;[a2(a3+a4)  =  -977,300 

14  =  ala2a3a4  1-  11.31 
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APPENDIX  M 


General  Equations  for  the  Pore  Air  Effect 

The  soil  phase  relationships  are  well  known.  Consider  the  element  of 
dry  soil  shown  below. 


The  void  ratio,  e,  is  the  ratio  of  void  volume  to  solid  volume. 


(M.l) 


whereas  the  porosity  is  the  ratio  of  void  volume  to  total  volume. 


n 


(M.2) 


Since  total  volume  equals  solid  volume  plus  void  volume, 

VT  =  V$  +  Vv  (M.3) 

Equations  (M.l)  and  (M.2)  can  be  written  in  the  forms 
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nV, 


e  = 


n  = 


Vs  =  VT  -  Vv  =  VT  -  nVT  =  rT_" 


eV 


(M.4) 


VT  '  Vs  +  Vv  ~  vs  +  eVs  " 


(M.5) 


For  one  dimensional  wave  propagation  and  pore  air  diffusion,  we  have 
the  following  situation: 


<r  •  / 


POK.E  A  IE 
PHASE. 


SO/t  S-KELETOhJ 
PHASE 


6 


* 


i 


& 's 


(r+iS-4-h  f 


JX  r/  ’  '  I  '  px 

The  volume  of  solids  in  a  given  soil  element  remains  constant,  so  that 

[Taylor  (1948:227)] 

VT0  -  dx-1  .  Vs  +  Vvo  =  Vs  *  e0Vs  =  »s(l  ♦  e0> 


(H.6) 


and  therefore. 


v  -  T0  _  dx’l 

vs  1  +  e_  1  +  e. 

o  o 


(M.7) 


The  equation  of  motion  for  the  soil  skeleton  is 


H  a*-1  -  f  "x-1  -  »svs  ^ 
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E 


where,  by  definition 


o  =  a  +  p 


Effective  stress  is  a  function  of  strain. 


=  o(e) 


strain  is  related  to  void  ratio. 


V-m  -  VT  (1  +  ejV„  -  (1  +  e)V 


and  void  ratio  is  therefore  related  to  strain. 


e  =  e  +  (1  +  e  )  — 
c  co  v  o'  ax 


(M. 13) 


For  adiabatic  flow. 


p  =  aP: 


(M. 14) 


The  equation  of  pore  air  flow  is 


1 1R 


(M. 15) 


and  the  equation  of  conservation  of  air  mass  is 


^  £».«(»- ffJf*'1 -- -it  W- 


I 
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or 


IT  ■ 


<»ae> 


(M.16) 


Inserting  Equation  (M.9)  in  Equation  (M.8)  yields 


2  1  +  e  - 

_3 _ U  ^  _ _ 0  /  3g  ^  8p  \ 

at2  =  '  »s  'ax  ax' 


(M.17) 


and  from  Equations  (M.10)  and  (M.ll)  we  obtain 


3 a  _  j3a  _3e  _  F  2 
3X  3e  3X  tT 

ax^ 


(M. 18) 


where 


p  3a 
tT  “  3c 


(M. 19) 


Thus,  Equation  (M.17)  can  be  written  in  the  form 


a2„  1  *  e 

3  U  l 


~  ^et  77  “  f  3 


(M.20) 


Inserting  Equation  (M. 15)  into  Equation  (M.16)  yields 


tVt!  (pae)  =  B1  if  ^alx* 


(M.21) 


Inserting  Equation  (M.14)  into  Equation  (M.21)  yields 
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or 


1  [>  3  (r^h)  4-  rjl/lf  3el 

l“e7  Le  IF  (p  >  p  Tt J  “ 


B,  2  -  +  1 

1  3  (PY  ) 


-  +  1  3X 
Y 


B,  2  —  +  1 

1  3  (PY  ) 


-+  1  ax' 

Y 


or 


3  <»1/T> 


1 


(n  +|H) 
'  0  3X ' 


B,  .2  -  +  1 

X  A  (PY 


(i  +  1)  ax2 
l'y  ' 


)  -  P 


1/y  3  u 
axat 


(M .22) 


An  alternate  form  of  Equation  (M.22)  is 


1 


V 

Y 


i  _  1 


i£ 

at 


1 


(n  +  — ) 
v  o  ax' 


B1  (pWY  7^  ‘  p 


1/y  32u1 

axat 


1 


.  au 
n  +  — 
o  ax 


V 

Y  M 


-  -  1 


(f): 


ax  J 


-  P 


1/y  a  u 
axat 


or 


3E  _ 

at 


1 


n  +  1H 
o  ax 


r 

o 

B1 

(f  >2  * 

a2u 

YP  axat  1 

k 

J 

(M.23) 


Equations  (M.20)  and  (M.22)  are  coupled,  nonlinear,  second  order 
partial  differential  equations  for  the  soil  skeleton  displacement,  u,  and 
the  pore  air  pressure,  p.  Gravitational  acceleration  does  not  appear 
explicitly  in  Equation  (M.20),  because  the  stresses  involved  are  assumed 
to  be  "live"  stresses,  i.e.,  stresses  in  excess  of  those  existing  under 
geostatic  conditions.  It  is  the  live  stresses  which  cause  motion.  Note, 
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however,  that  when  live  stresses  are  used,  they  must  be  expressed  in  terms 
of  live  strains.  This  means  that  if  the  stress-strain  relation.  Equation 
(M.10),  is  nonlinear  the  geostatic  stress-strain  condition  must  be 
determined,  because  it  must  be  used  as  the  point  of  departure  for  the  live 
stress-live  strain  relation. 

The  effective  permeability,  B^,  has  been  assumed  constant  in 
Equations  (M.15),  (M.21),  (M.22)  and  (M.23),  although  is  known  to 
vary  with  strain.'  This  assumption  appears  reasonable  prior  to  spall, 
because  the  strain  required  to  cause  particle  separation  is  small,  and 
because  the  exact  relation  between  effective  permeability  and  strain  is 
often  not  precise. 
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Table  3.1.  One-Dimensional  Calculations 


Calc. 

No. 

Description 

Tensile  Cutoff 

P  •  ,  MPa 
min 

Peak  Incident 
Velocity,  m/s 

Comments 

1 

no  spall  allowed 

100. 

1.0 

3  cycles  to 
spall 

2 

baseline  spall 
calculation 

0.10 

1.0 

3  cycles  to 
spall 

3 

load  variation 

0.10 

2.0 

3  cycles  to 
spall 

4 

no  tensile 
strength 

0.00 

1.0 

3  cycles  to 
spall 

5 

higher  tensile 
strength 

0.30 

1.0 

3  cycles  to 
spal  1 

6 

immediate  spall 

0.10 

1.0 

1  cycle  to 
spall 

7 

forced  rejoin 

0.10 

2.0 

1.0  MPa 

O.P.  with 
expon. 
decay  at 

300  ms 
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SUMMARY  OF  PEAK  PORE  AIR  PRESSURE  BELOW  THE  SURFACE  OF  SAND 
DURING  EXPOSURE  TO  A  100  PSI  SHOCK  IN  AN  8-INCH  SHOCK  TUBE 


t 


<3* 


co 

rv. 

<r» 


<D 


* 

o 


S- 

<V 


TABLE  4.2 


TABLE  4.3 


CALCULATED  PEAK  PORE  AIR  PRESSURE  AND  TIME 
OF  OCCURRENCE  FOR  ZERNOW'S  8- INCH  VENTED  SAMPLES* 


GROUP 

D 

e 

X 

e 

PMAX 

tmax 

tMAX 

CM2/SEC 

IN 

PS  I 

SEC 

I 

3.51X105 

0.463 

33.0 

0.465 

46.1 

0.35 

0.032 

II 

3.51X105 

0.463 

57.6 

0.811 

16.2 

0.40 

0.037 

III 

3.51X105 

0.639 

60.0 

0.845 

12.3 

0.37 

0.034 

IVa 

5.84X105 

0.279 

33.0 

0.465 

48.1 

0.39 

0.022 

IVb 

5.84X105 

0.279 

57.6 

0.811 

16.9 

0.43 

0.024 

*Note: 


r  x  x 

?  -  T"  7T 

l2TMflx  <180.34)2Thax 
tMAX  D  D 
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TABLE  4.4 


ADJUSTED  LINEAR  DIFFUSION  PARAMETERS  FOR 


ZERNOW'S 

8-INCH  SHOCK  TUBE  VENTED 

SAMPLES* 

GROUPS 

Bi 

n 

D 

td 

a 

3 

CM4/DYNE  SEC 

cm2/sec 

SEC-1 

I, II 

3.4X10-3 

1/3 

7.03X104 

2.16 

5.0 

2.31 

III 

3.4X10-3 

1/3 

7.03X104 

2.16 

6.9 

3.19 

IV 

6.0X10"3 

0.85/2.4 

1.17X105 

3.59 

5.0 

1.39 

*Note: 


B.P  6.89  X  106  B 


=  °^D  =  D(l.O) 

'  l2  (180. 34) 2 

l2ot  _  ( 180. 34)Za 
D  D 
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TABLE  4.5 


CALCULATED  PEAK  PORE  AIR  PRESSURE  AND  TIME 
OF  OCCURRENCE  FOR  ZERNOW'S  8-INCH  VENTED  SAMPLES, 
USING  ADJUSTED  LINEAR  DIFFUSION  PARAMETERS* 


GROUP 

D 

B 

X 

e 

PMAX 

tmax 

tMAX 

CM2/SEC 

IN 

PSI 

SEC 

I 

7.03X104 

2.31 

33.0 

0.465 

34.4 

0.21 

0.097 

II 

7.03X104 

2.31 

57.6 

0.811 

12.0 

0.26 

0.120 

III 

7.03X104 

3.19 

60.0 

0.845 

8.9 

0.24 

0.111 

IVa 

1.17X105 

1.39 

33.0 

0.465 

38.9 

0.25 

0.069 

IVb 

1.17X105 

1.39 

57.6 

0.811 

13.6 

0.30 

0.083 

*Note : 


r  _  x  _  x 
K  1  "  71 


tMAX 


]2tmax  _  <180-34>2tmax 
D  5 
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TABLE  4.6 


ADJUSTED  LINEAR  DIFFUSION  PARAMETERS  FOR 
8-INCH  SHOCK  TUBE  UNVENTED  SAMPLES* 


GROUPS 

Bi 

n 

D 

To 

a 

B 

cm4/dyne  sec 

cm2/sec 

SEC"1 

1,11 

3.4X10'3 

1/3 

7.03X104 

0.540 

5.0 

9.25 

III 

3.4X10'3 

1/3 

7.03X104 

0.540 

6.9 

12.77 

IV 

6.0X10'3 

0.85/2.4 

1.17X105 

0.897 

5.0 

5.56 

*Note : 


B, P  6.89X106  B 


=  Dt°  =  D(l.O) 

D  4T7  (360. 68)2 

_  41 2a  _  (360.68)2a 
D  D 
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TABLE  4.7 


CALCULATED  PEAK  PORE  AIR  PRESSURE  AND  TIME 
OF  OCCURRENCE  FOR  8-INCH  UNVENTED  SAMPLES, 


USING 

ADJUSTED 

LINEAR 

DIFFUSION 

PARAMETERS* 

GROUP 

D 

8 

X 

i 

PMAX 

tmax 

*MAX 

CM2/SEC 

IN 

PS  I 

SEC 

I 

7.03X104 

9.25 

33.0 

0.232 

70.9 

0.010 

0.019 

II 

7.03X104 

9.25 

57.6 

0.406 

44.2 

0.050 

0.093 

III 

7.03X104 

12.77 

60.0 

0.423 

38.8 

0.045 

0.083 

lVa 

1.17X105 

5.56 

33.0 

0.232 

73.4 

0.010 

0.011 

IVb 

1.17X105 

5.56 

57.6 

0.406 

51.5 

0.085 

0.095 

*Note: 


r  _  x  _  x 

K  2T  W 

_412Thax  (360.68)2Tmax 
xMAX  D  D 
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TABLE  5.1 


Scaled  Airblast  Parameters,  Large  Yield 
Pacific  Nuclear  Data  (lOkt  -  lOmt) 


Scaled  Range 

"R,  km/(kt)l/3 

Peak  Positive 
Overpressure 
aP+,  kPa 

Peak  Negative 
Overpressure 
aP~,  kPa 

0.134 

345 

27.6* 

0.151 

297 

15.2 

0.194 

159 

13.8 

0.244 

94.5 

12.4 

0.933 

10.0 

1.86 

0.956 

10.1 

0.90 

0.977 

10.3 

2.74 

2.948 

2.21 

0.66 

2.996 

1.65 

0.56 

4.060 

1.17 

0.44 

0.048 

1131 

101? 

0.207 

175 

21.9 

0.244 

130 

13.8 

0.394 

48.3 

9.03 

0.580 

24.4 

4.34 

2.063 

2.62 

0.90 

2.063 

3.10 

1.17 

0.205 

121 

14.4 

0.281 

77.2 

11.6 

0.381 

60.8 

9.10 

0.221 

128 

14.8 

0.298 

85.2 

8.62 

0.662 

18.7 

4.00 

1.042 

9.30 

2.17 

1.590 

3.83 

0.97 
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EXPERIMENTAL  SPALL  DATA  PROM  NEAR-SURFACE  DETONATIONS 
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TABLE  6.2 

Maximum  Scaled  Radius  of  Spall, ~RS, 
for  Various  Explosive  Charge  Configurations 


• 

Charge 

Configuration 

Rs  (m/100  ton1/3) 
with  Negative 
Airblast  Wing 

Rs  (m/100  ton1/3) 
without  Negative 
Airblast  Wing 

• 

TNT 

HOB 

No  Data 

32 

TNT 

STS 

125-200 

40-  80 

TNT 

STC  or  HBS 

140-200 

90-105 

TNT 

Berm 

N/A 

140 

<t 

Nuclear 

DOB 

N/A 

48 

ft 


'ft 


'ft 
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Depth  in  Meters 


Material  Models 


Zoninc 


0  fO  mf 


Free  Surface 


Soil 

p  =  2.0  g/cc 
u  =  0.30 


Rock 

p  =  2.5  g/cc 
u  =  0.20 


60  zones 
at  0.1  m 


60  zones 
at  0.15  m 


60  zones 
at  0„3  m 


Soi  1 : 


C,  =  561 

L  V  m/s. 


■M  = 

Cu  =  1122  m/s 


°MIN»  vaHed 


Axial  Strain, 


M2  =  M  = 

C.  =  Cu  =  2960  m/s 


Axial  Strain,  % 


Transmitting  Boundary  with 
Incident  Velocity,  V(t): 


Linear  Rise 


Sinusoidal  Decay 
.  Zero, Incident  Velocity 


Time  (ms) 


Fig.  3.2.  One  Dimensional  (Uniaxial  Strain)  Spall  Calculations 


6 


Velocity  (m/s) 


Fig.  3.3.  Application  of  Gravity  to  Grid  with  Dynamic  Relaxation 


Stress  (MPa) 


Fig.  3.4.  In  Situ  (Gravity)  Stress  and  Pressure  Distribution 


New  Time  =  tn 
Old  Time  = 

Current  Time  =  (tn  +  t^)/2 
Timestep  =  At  =  tn  -  t^ 

Note:  Assuming  small  strains 


Fig.  3.5.  Transmitting  Boundary  with  Incident  Velocity 


Depth  (m)  Depth  (m) 


Time  (ms) 


a.  Characteristic  Plane  with 
Pressure  Boundary  at  27  m 
(No  Tensile  Failure) 


Time  (ms) 


b.  Characteristic'  Plane  with 
Transmitting  Boundary  at  27  m 
(No  Tensile  Failure) 


Fig.  3.6.  Characteristic  Planes  for  Two  Boundary  Conditions 


0.20  r  i  •  -  =  Compression 


Calc.  #1 
No  Spall 


Fig.  3.8.  Spall  and  No  Spall  Time  History  Comparisons  at  Depth  =  0.5 


Soil 


0.50 


Position  Time 


FIGURE  4.1a 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  6  =  0.463,  T  =  0.05 


Pressure,  PSI 


01 


FIGURE  4.1b 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  B  =  0.463,  T  =  0.15 


Pressure,  PSI 


Normalized  Position 


FIGURE  4.1c 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  g  =  0.463,  T  =  0.25 
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Normalized  Position 


FIGURE  4.2a 

Dynamic  Pore  Pressure  Isochrone  for 
Vented  Sample;  6  =  0.639,  T  =  0.05 


ssure,  PSI 


! 

i 

I 


Normalized  Position 


FIGURE  4.2b 

•  Dynamic  Pore  Pressure  Isochrone  for  a 

Vented  Sample;  B  =  0.639,  T  =  0.15 


175 


Normalized  Position 


FIGURE  4.2c 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  8  =  0.639,  T  =  0.25 
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80 


<• 


FIGURE  4.3a 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  8  =  0.279,  T  =  0.05 


Pressure,  PSI 


FIGURE  4.3b 

•  Dynamic  Pore  Pressure  Isochrone  for  a 

Vented  Sample;  8  -  0.279,  T  =  0.15 


1.25 


178 


*• 


FIGURE  4.3c 

Dynamic  Pore  Pressure  'sochrone  for  a 
Vented  Sample;  3  =  0.279,  T  =  0.25 
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Normalized  Time 


FIGURE  4.4a 

Dynamic  Pore  Pressure  Response  in  a 
Vented  Sample;  6  =  0.463,  £  =  0.465 


17.5 


Normalized  Time 


FIGURE  4.4b 

Dynamic  Pore  Pressure  Response  in  a 
Vented  Sample;  3  =  0.463,  E,  =  0.811 
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FIGURE  4.5 

Dynamic  Pore  Pressure  Response  in  a 
■Vented  Sample;  0  =  0.639,  £  =  0.845 
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microcopy  resolution  test  chart 

NATIONAL  BUREAU  OF  STANDARDS-  1963-A 
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Normalized  Time 


FIGURE  4.6a 

Dynamic  Pore  Pressure  Response 
Vented  Sample;  6  =  0.279,  £  =  0 


in  a 
.465 


Normalized  Time 


FIGURE  4.6b 

Dynamic  Pore  Pressure  Response  in  a 
Vented  Sample;  6  =  0.279,  £  =  0.811 


0  0.25  0.50  0.75  1.00 

Normalized  Position 


FIGURE  4.7a 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  g  =  2.31,  T  =  0.05 
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1.25 


85 


0.00  0.25  0.50  0.75  1.00 

Normalized  Position 


FIGURE  4.7b 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  6  =  2.31;  T  =  0.15 


<% 


1.25 
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Normalized  Position 


FIGURE  4.7c 

Dynamic  Pore  Pressure  Isochrone  for 
Vented  Sample;  B  =  2.31;  T  =  0.25 


Normalized  Position 


FIGURE  4.8a 

Dynamic  Pore  Pressure  Isochrone  for 
Vented  Sample;  6  =  3.19;  T  =  0.05 


Pressure,  PSI 


Normalized  Position 

FIGURE  4.8b 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  8  -  3.19;  T  =  0.15 
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Pressure,  PSI 


FIGURE  4.8c 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  B  =  3.19;  T  =  0.25 


1.25 
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Normalized  Position 


FIGURE  4.8d 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  6  =  3.19;  T  =  0.35 


Normalized  Position 


FIGURE  4.8e 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  B  =  3.19;  T  =  0.45 


FIGURE  4.8f 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  B  =  3.19;  T  =  0.55 
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Normalized  Position 


FIGURE  4.9a 

Dynamic  Pore  Pressure  Isochrone  for 
Vented  Sample;  8  =  1.39;  T  =  0.05 


Normalized  Position 


FIGURE  4.9c 

Dynamic  Pore  Pressure  Isochrone  for  a 
Vented  Sample;  B  =  1.39;  T  =  0.25 
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FIGURE  4.10a 

Dynamic  Pore  Pressure  Response  in  a 
Vented  Sample;  $  =  2.31;  £  =  0.465 
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Normalized  Time 


FIGURE  4.10b 

Dynamic  Pore  Pressure  Response  in  a 
Vented  Sample;  0  =  2.31;  £  =  0.811 
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« 


FIGURE  4.11 

Dynamic  Pore  Pressure  Response  in  a 
Vented  Sample;  6  =  3.19,  £  =  0.845 
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FIGURE  4.12b 

Dynamic  Pore  Pressure  Response  in  a 
Vented  Sample;  3  =  1.39;  £  =  0.811 
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Pressure,  PSI 


FIGURE  4.13c 

Dynamic  Pore  Pressure  Isochrone  for  an 
Unvented  Sample;  g  =  9.25,  T  =  0.15 


Normalized  Position 


FIGURE  4.14b 

Dynamic  Pore  Pressure  Isochrone  for  an 
Unvented  Sample;  0  =  12.77,  T  =  0.10 
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FIGURE  4.14c 

Dynamic  Pore  Pressure  Isochrone  for  an 
Unvented  Sample;  g  =  12.77,  T  =  0.15 
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Pressure,  PS 


FIGURE  4.15a 

Dynamic  Pore  Pressure  Isochrone  for  an 
Unvented  Sample;  0  =  5.56,  T  =  0.05 


Dynamic  Pore  Pressure  Isochrone  for  an 
Unvented  Sample;  6  =  5.56,  T  =  0.15 
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FIGURE  4.15c 

Dynamic  Pore  Pressure  Isochrone  for  an 
Unvented  Sample;  B  =  5.56,  T  =  0.25 
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Figure  4.16a 

Dynamic  Pore  Pressure  Response  in  an 
Unvented  Sample;  6  =  9.25,  £;  =  0.232 


Normalized  Time 


FIGURE  4.16b 

Dynamic  Pore  Pressure  Response  in  an 
Unvented  Sample;  6  =  9.25,  £  =  0.406 


Dynamic  Pore  Pressure  Response  in  an 
Unvented  Sample;  6  =  12.77,  £  =  0.42 


Pressur 


FIGURE  4.18a 

Dynamic  Pore  Pressure  Response  in  an 
Unvented  Sample;  8  =  5.56,  £  =  0.232 
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0.01  0.10  1.0  10 

Scaled  Range,  R(km/kt1/3)nuclear 

FIGURE  5.2 

Peak  Negative  Airblast  Overpressure  Versus  Scaled  Range,  Nuclear 


x  Negative  Airblast 
Wing  Region 


/ 


z  =  depth  of  spall  at  specified 
range 

z$  =  maximum  depth  of  spall 

R  =  maximum  radius  of  spall  for 
a  depth  of  0.5  m  (arbitrary) 


FIGURE  6.1 

Sketch  of  Typical  Spall  Regions 
and  Nomenclature 
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FIGURE  6.2 

Estimated  Extent  of  Spall,  PHG  1-07 


FIGURE  6.3 

Typical  Vertical  Velocity  Rejoin  Amplitude,  PRE-MINE  THROW  IV 


Soil  Tensile  Strength,  oj  (kPa) 


Scaled  Range,  R(m/0.5  ton 
FIGURE  6.4 

Summary  of  Vertical  Velocity  Rejoin  Amplitude  Versus  Scaled  Range,  PRE-MINE  THROW  IV 
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Maximum  Spall  Depth,  zg  (m) 

FIGURE  6.6 

Maximum  Depth  of  Spall  Versus  Depth  to  the  Second  Layer 
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Yield,  W  (tons) 
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zs  =  3.2  (W)1/6,m 


Maximum  Spall  Depth,  z$  (m) 
FIGURE  6.7 

Maximum  Depth  of  Spall  Versus  Yield 
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Scaled  Range,  R  (km/kt  '  ) 

FIGURE  7.1 

Spall  Depth  in  the  Negative  Airblast  Wing  Region 


